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/*MDL_BEGIN 

MDL_SOURCE_FILE_NAME : RI CKS_R9 6 MATH . H 

PKG_ACRONYM : R96MAI MDL_SOURCE_VERSION# : VI . 0 

MDL_DESCRIPTION : 

This module is an include file that provides some of the 

math support for TTI's REAL 9 6 math datatype. 
******************************************************************* 

* Copyright (c) 2001 by * 

* Touch Technologies Inc., San Diego, California. * 

* ALL RIGHTS RESERVED. * 
**************************************************************************** 

R96 datatype example implementation 
Data stored in Little endian 

bits 0-63 are the Integer Part (IP) (signed, so 63 bits of precision) 

bits 64-95 are the Fractional Part (FP) (signed, modulus 100,000,000 - ie. 8 digits) 



Sign is generally stored in both IP and FP, however, 
determine the sign of a r96, as follows: 

If IP != 0 then IP is 2s complement and can be used to determine the sign 
If IP == 0 then FP is 2s complement 

Because zero has no sign, to determine the sign if IP=0, you have to 
look at the sign stored in fractional part. 

/* In specific embodiments, routines are provided for unsigned mutliprecision 



Tti 

03 



Id 

pas 
III 



Integer Math (64, 96, and 128 bit) . _ 

MACROS and routines are organized herein so that 64, 96, and 128 bit variants are in order for each 
routine or MACRO . This makes maintenance easier because the variants are very similar. 

Most math is done in unsigned 96bit integer format (scaled). Conversion to 96bit unsigned integer state is 
done by removing the sign and then multiplying the IP by 100,000,000 and adding in the FP. For multiply 
and divide, some 128 bit math may be required to prevent bit loss. Optimizations have been performed to 
minimize the number of bits required, especially when multiplication is required. 

There are MACROS to perform 196 compares, additions and subtractions. Most MACROS assume that the 

parameters are of type (unsigned long *) _ . 

Miscellaneous MACROS 



iviisceiian epin> MMv*rwa> _ _ . 

There are MACROS to perform comparisons such as «, 1 =,>,<, and compares against ZERO. There are also 
MACROS to copy an integer and to negate an integer. . 



Bit shift MACROS 



There are MACROS to perform right and left bit shift with single bit shifts and counted shifts. Counted 
shifts MUST be 31 bits or less. There are MACROS to check for overflow as we ll used in multiplication. — 

Addition MACROS: 



If a 2 operand ADD is needed then the SUM and the second ADDEND must be the same variable and NOT the 
first ADDEND, i.e. a - a + b is coded as add (b, a, a) . 

There are MACROS to handle overflow as well necessary for multiply. The i##add32 variant adds a 32 bit 

integer to the appropriate extended precision integer providing a shortcu t for that case 

Subtraction MACROs: - 



These MACROS perform unsigned subtraction. 



Multiplication Routines 



The basic mechanism is to determine whether the multiplicand or the 
multiplier is smaller and use that number as the multiplier. The 
multiplier is then shifted left until bit ZERO is set. Each time 
the multiplier is shifted left, the multiplicand is shifted right. 
When bit ZERO is set in the multiplier, the current shifted 
multiplicand is added to the product. When the multiplier becomes 
ZERO, the multiplication stops and the product is now complete. 

We do the multiplication by BIG "switch" statement that does 8 

bits at a time to save time over doing a single bit at a time. 

The "ixxmultiply" and "ixxmultiply_oiverf low" routines are created by 

a program named "create_ixxmultiply_r96math.c" . It does the 64, 96, 

and 128 bit routines. 

There are also routines to test for and handle overflow. Overflow is 
identical to non-overflow except the shift and add routines called 
check for the overflow condition. 

There are also some special MACROS that multiply specifically by 10, 
10,000, and 100,000,000. These speed up the various R96 conversions 
and operations. 
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i NOTE: These routines are 3 operand routines!!! One CANNOT use the same variable for the ! 

product as either the multiplicand or the multiplier!!! 

1 

! There is a set of routines to locate the highest set bit. These are 

! used to optimize overflow handling in multiplication. It turns out that 

! overflow can be easily predicted by adding the highest set bits and 

! comparing them against the number of the sign bit. >= will be overflow. 

The final result is also compared to see if the sign bit is set which 
implies overflow before sign handling is done. 



1 Integer Divides 



The algorithm shifts the divisor over until it is >= the dividend and 
then shifts it back one bit. The divisor is then subtracted from the 
dividend. The quotient gets the value of the number of bit shifts that 
were done added to it. The difference of the subtract becomes the 
remainder. This process is repeated until the remainder is less than 
the divisor. 



R96 miscellaneous routines and MACROS 



r96abs is a MACRO that converts an R96 number into a positive value 
r96negate is a MACRO that negates an R96 number 

r96copy is a MACRO that copies an R96 number from one location to another 
r96cvtr96toi96 is a routine that converts an R96 number into an unsigned 
96 bit integer. The sign is returned as the function value (0 if positive 
and 1 if negative). The IP is multiplied by 100,000,000 and then the 
FP is added in. 

r96cvtr96toi96_4digitfp is a special case routine of r96cvtr96toi96 that 
converts an R96 number into an unsigned 96 bit integer when modulus 10 , 000 
of the FP == ZERO. The sign is returned as the function value {0 if 
positive and 1 if negative). The IP is multiplied by 10,000 and then the 
FP is added in. 

r96cvti96tor96 is a routine that converts an unsigned 96 bit integer 
into an R96 number. The sign is a parameter. The 196 value is 
divided by 100,000,000. The quotient becomes the IP and the remainder is 
the FP. 

r96cvtasctobin converts an ASCIZ string into an R96 number. 

r96cvtbintoasc converts an R96 number into an ASCIZ string. 

If this routine is to be used as part of the TTI R96 math stuff, 

it should have an optimization that uses a table lookup where the 

the binary is divided by 10,000 converting 4 digits at a time, 

instead of 10 which only converts 1 digit at a time. 



R96 division 



Divide is accomplished by converting the dividend and divisor to 96bit 
integers and then doing a 96bit divide. The remainder is then converted 
to a 128 bit integer and multiplied by 100,000,000. This value is then 
divided by the original divisor using the least required bits to get the 
FP portion of the result quotient. Overflow is tested by checking bit 
63 of the result quotient. Overflow can happen when dividing by a 
1 fractional entity, i.e. 100 / .1 == 100 * 10. 
i 

i we call r96divide_nohilong as an optimization when the IP part of both 
! numbers is < 2*32. This routine is able to save 32 bits of precision 
I in the operations. 

Overflow handling is done by checking the unsigned final quotient 
(IP portion) against 0x7FFFFFFFFFFFFFFF before signs are applied. 
If the value is larger, then runtime$error is called to signal the 
error. The runtime$error in the R96 math .H files is intended 
to be over- written. The definition is conditional. 



R96 multiplication 



The multiplication starts with a conversion to 128 bit unsigned 
integers. a 128 bit multiply is then done. The product is 
then divided by (100,000,000 * 100,000,000) to get the IP part 
of the product. The remainder of that divide is divided by 
100,000,000 to get the FP part of the product. 

We start by checking for possible optimizations. If both FPs == ZERO 
then we call r96multiply_nofp which performs 164 math. If either 
the multiplicand or the multiplier == ZERO, then we call 
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i r96multiply_nofp_multiplicand or r96multiply_nofp_multiplier as 

! appropriate. Finally, if modulus 10,000 of both FPs == 0 then we call 

! r96multiply_4digitfp. The last 3 optimization routines perform 196 oit 

! math instead of 1128 bit thus running faster. 

j 

i Overflow handling is done by checking each left shift and addition for 
! overflow as well as checking the unsigned final product (IP portion) 
! against 0x7FFFFFFFFFFFFFFF before signs are applied. If the value is 
! larger, then runtime$error is called to signal the error and processing 
I continues until a solution is arrived at. The runtime$error in the R96 
! math .H files is intended to be over-written. The definition is 
i conditional. THIS IMPLIES THAT runtime$error CAN BE CALLED MULTIPLE 
I TIMES, Noticing an overflow does not directly stop the processing of 
I the multiplication. 
*/ 

#include "ixxmath_nomultiply.h" 
#include "ixxmath_multiply.h" 
#include "r96math.h" 
/ *MDL_END*/ 



* COPYRIGHT * 

* Copyright (C) 2001 * 

* by Touch Technologies Inc., San Diego, California * 

* This software is furnished under a license and may be used and * 

* copied only in accordance with the terms of such license and with * 

* the inclusion of the above copyright notice. This software or any * 

y * other copies thereof may not be provided or otherwise made * 

|pj * available to any other person. No title or ownership of the * 

M * software is hereby transferred * 

********************************************************** 

.f| FACILITY: SheerPower Language 

? !| ABSTRACT : 

%4 This module provides the real integer 18.8 math class for 

yrf both Alpha & Intel Pentium Processors.*/ 

Module for the real integer 18.8 math class for Alpha & Intel Pentium 

#include "c-universal .h" 
#include <stdio.h> 
ttinclude <string.h> 
#include "c-intmsg.h" 
#include <time.h> 
# include <math.h> 
#include "ricks_r96math. h" 
IH- typedef struct 

O { int wholedigits; /* number of whole digits eg. 99.678 has 2 whole digits*/ 

HI int fractdigits; /* number of fractional digits eg. 99.678000 has 3 fractional digits*/ 

int neg; /* if the number is negative*/ 

int rlen; /* result length eg. 99.678 has an rlen of 6*/ 

} real info; 
long ran_random ( ) ; 
int randomize () ; 
void scaledown (void *_a, int n) ; 
int ator{real *a, char *b> ; 

void rtoa (real *, void *, int f digits, realinfo *rinfo, int flags); 
inline void scaleup (void *_a) ; 
static real scale = 1; 

#define LARGEST NUMBER (~ ( ( (int64) 1) «63) ) ,.,,,,,,,,,,,,,,,,,,/ 
7 /////////-/ / i t // / // / / / t f / / ///// t / t ////// / /// / t 11 f ' f 1 1 f ' f 1 1 ' 1 { 1 ' 1 ' 

n MTHSDEXP DOUBLE EXPONENTIAL ______ 

II Brief description: Return the exponential of the input value. 
// Expected: d = the address of a 'double' 

// Result: return = the exponential of 'd' as a double ,,,,,,,,,, 

double routine mth$dexp (double *d) 

{ return exp (*d) ; // call exp() passing the value of «d'} ,,,,,,,,,, 

} / 7 / 77 /7 / /////////////// / ////////////// > > > > > > t—L ' i ' ' t i ' l ' ' 1 > 1 > ' 1 1 ' ' ' 

I /; MTH$DLOG2 DOUBLE BASE 2 LOGARITHM 

I : : : ; : : : ! : : ; ; : 7 5 7 7 ; 7 7 7 7 T~ 



//,//////////////////////////////////////// ri / t f i > i / i i ii / / / / / / / / * * > 
II Brief description: Return the Base 2 logarithm of the input value. 
// Expected: d = the address of a 'double' 

// Result: return = the base 2 logarithm of 'd' as a double 
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///////////////////////////////////////////////>/////////////////// 
double routine mth$dlog2 (double *a) 
{ return log (*a) * 1.442695099;} 

///////////////////////////////////////////////'/////////////////// 



„ MTH$DN1NT DOUBLE TO INTEGRAL PART 

/////////////////////////////////////////////////////////////////// 
// Brief description: Remove the fractional part from a double and return result. 

// Method: Convert the double to a large 64-bit integral value. This removes 

// the fractional part. Then convert the large 64-bit integer back 

// to a double, and return that value. 

// Expected: a = the address of a 'double' 

// Result: return * the integral part of 'a' as a double. 

/////////////////////////////////////////////////////////////////// 
double routine mth$dnint (double *a) 

{ int64 i = (int64) * a; // convert to a large interger 

double d = (double) i; // then to a double 

return d; // return integral part as a double} 

//////////////////////////////////////////////// /////////////////// 



u OTS$POWDD RAISE A DOUBLE TO A DOUBLE EXPONENT 

///////////////////////////////////////////////////////////// ////// 
// Brief description: Raise a double 'x' to a double exponent 'y'. 



D 



^5? 



// Expected 
II 

f I Result 



x = the base double 

y = the exponent double 

return = double 'x' raised by 'y' 

/////////////////////////////////////////////// / / // // // // // // // // // 
double routine ots$powdd (double x, double y) 
{ double r = pow(x, y) ; 
return r; } 

/////////////////////////////////////////////////////////////////// 

f OTS$POWJJ RAISE AN INTEGER TO AN INTEGER EXPONENT 

///////////////////////////////////////////////,//////////// /////// 
II Brief description: Raise an integer 'x' to a integer exponent 'y' . 
// Expected: x = the base integer 

// : y = the exponent integer 

// Result: return - double 'x' raised by 'y' 

v f\ /////////// / / / / / / / / / // // / / / // // // // // / / / // // // // // // // // // // / / / // / / 
*JJ[ double routine ots$powjj (int _x, int _y) 
{ double y = _y; 

if! double x = _x; // convert to double's first for pow() 

2^ double r = pow(x, y) ; 

s ¥ return r ; } 

/////////////////////////////////////////////////////////////////// 
a OTSSPOWDJ RAISE A DOUBLE TO AN INTEGER EXPONENT 

3 '///////////////////////////////////////////////< / // // // // // // // // // 

// Brief description: Raise a double 'x' to an integer exponent 'y' . 



if Expected 

W // 

// Result 



x = the base double 

y = the exponent integer 

return = double 'x' raised by 'y' 
/////////////////////////////////////////////////////////////////// 
p| double routine ots$powdj (double x, int _y) 
JsJ { double y = _y; // convert to double first 

double r = pow(x, y) ; 
ilj return r;} 

////////////////////////////////////////////// / t ////////////////// / 

h MTHSDSQRT SQUARE ROOT OF A DOUBLE 

/ / i /////////////////////// t / / / / '/ /"/""/ ////// / / / / / / t ////////////////// / 

// Brief description: The square root of a double. 

// Expected: d = the address of the double 

// Result: return = the square root of ! d' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dsqrt (double *d) 
{ return sqrt { *d) ; } 

///////////////////////////////////////////////;/////////////////// 

n MTHSDLOG NATURAL LOGARITHM 

/////////////////////////////////////////////////////////////////// 

// Brief description: The natural logarithm of a double as a double. 

// Expected: d = the address of the double 

// Result: return = the natural logaithm of 'd' as a double 

t ////////////////////////// i /////////////////// t //////////// i ii / / / / 
double routine mth$dlog (double *d) 
{ return log (*d);} 

t f l l i f I t I I t ( t t f I I I t t f ! I I ! I I l f t l I l l 1 I l i I t J I I ! i I f > t I I t t t t t I ' ' ' I t I l l ' I 

h MTHSDLOG10 COMMON LOGARITHM 

////////////////////////////////////////// /////>/////////////////// 

// Brief description: The common logarithm of a double. 

// Expected: d = the address of the double 

// Result: return = the common logarithm of 'd' as a double 

///////////////////////////////////////////////*/////////////////// 
double routine mth$dlogl0 (double *d) 
{ return loglO (*d) ;} 
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/// / / / // // // // // // // // i /////// i /////////////////////////////////// i 



MTH$DSIN THE SINE OF ITS RADIAN ARGUMENT 



i f t i i t i i i f / / i / / / / // / // / / // / / // // // //'///////////// ////////////////// 

// Brief description: return the sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = the return the sine of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dsin (double *d) 

{ return sin (*d) ; // return the sine of its radian argument} 

i t i i i ' i i i i i t l I I l I l I I I I i I I I I I I I I { I I I f I I I I I t i I I l I I I ! I J I f i i I J I I I I i i i t f 



MTHSDSINH THE HYPERBOLIC SINE OF ITS RADIAN ARGUMENT 

i ////// f / / / / i iii / / / / // // // // / / / // // / / / / // // // / / / // // // // // // // // // / 

// Brief description: return the hyperbolic sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the hyperbolic sine of its radian 

// argument 'd 1 as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dsinh {double *d) 

{ return sinh (*d) ; // return the hyperbolic sine of its radian argument} 

/ / // // / ////// / /////////////////////// t ///////////(//// t / / / / // // / / / / 



// MTHSDCOS THE COSINE OF ITS RADIAN ARGUMENT ^ 

///////////////////////////////////////////////// ////////////////// 

// Brief description: return the cosine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the cosine of its radian 

// argument 'd 1 as a double 

t f t / f /////// f ///// f // f f /////// f /// f ////////////{ J ///////!/////,/// , 

double routine mth$dcos (double *d) 

{ return cos (*d) ; // return the cosine of the radian argument} 

///////////////////// / // / //// / / / / / / / //// / /?// / / t //// / ///,,//,//,// / 



MTH$DCOSH THE HYPERBOLIC COSINE OF ITS RADIAN ARGUMENT 



////////////////////////////////////////////// / < / // // // // // // // // / / 

// Brief description: return the hyperbolic cosine of its radian argument: 

// Expected: d = the address of the double 

// Result: return = return the hyperbolic cosine of its radian 

// argument ' d' as a double 

/ / // // / / / / // // // // // / ////////////////////////// * ////////////////// t 
double routine mth$dcosh (double *d) 

{ return cosh (*d) ; // return the hyperbolic cosine of its radian argument} 

///////////////////////////////////////////////>/////////////////// 



MTH$DTAN THE TANGENT OF ITS RADIAN ARGUMENT 

/////////////////////////////////////////////////////////////////// 

// Brief description: return the tangent of its radian argument 

// Expected: d = the address of the double 

// Result: return * return the tagent of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dtan (double *d) 

{ return tan (*d) ; // return the tangent of its radian argument} 

/////////////////////////////////////////////////////////////////// 
// MTHSDTANH THE HYPERBOLIC TANGENT OF ITS RADIAN ARGUMENT 

///////////////////////////////////////////////,-/////////////////// 
// Brief description: return the hyperbolic tangent of its radian argumert 
// Expected: d = the address of the double 

// Result: return = return the hyperbolic tangent of its radian 
// argument 'd' as a double 

/ / // / / / / / / / // // // / / // // // / / / // / // // // // // // // // i/ // // // // // // // // / / 
double routine mth$dtanh (double *d) 

{ return tanh (*d) ; // return the hyperbolic tangent of its radiar argument} 

////////// f / f ////// f // f //////////// f //////// t //,///// f // f // l I // f // l 

// MTHSDATAN THE ARC TANGENT OF ITS RADIAN ARGUMENT 

/////'///////////////////////////////////////,/,/////////////////// 

// Brief description: return the arc tangent of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the arc tagent of its radian 

// argument 'd 1 as a double 

////////// t // / /////////////////////////////// j //////////////////// / 
double routine mth$datan (double *d) 

{ return atan (*d) ; // return the arc tangent of its radian argument} 

/////////////////////////////////////////////////////////////////// 
n MTHSDASIN THE ARC SINE OF ITS RADIAN ARGUMENT 

/////////////////////////////////////////////////////////////////// 

// Brief description: return the arc sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the arc sine of its radian 

// argument 'd' as a double 

///////////////////////////////^////////////////y////////////////// 
double routine mth$dasin (double *d) 

{ return asin (*d) ; // return the arc sine of its radian argument } 

/////////////////////////////////////////////////////////////////// 
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n MTH$PACOS THE ARC COSINE OF ITS RADIAN ARGUMENT 

////////////// / / // // / / // // // // // // // / / /////////////// } >/ / J / // // // / / 

// Brief description: return the arc cosine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the arc cosine of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dacos (double *d) 

{ return acos (*d) ; // return the arc cosine of its radian argument} 

/////////////////////////////////////////////// f ////////////////// / 



K MTHSDABS THE ABSOLUTE VALUE OF A COMPLEX NUMBER 

/////////////////////////////////////////////// t ////////////////// / 

II Brief description: return the absolute value of a complex number 

// Expected: d = the address of the double 

// Result: return = return the absolute value of a complex 

// number as a double 

///////////////////////////////////////////////'/////////////////// 
double routine mth$dabs (double *d) 
{ return fabs (*d);} 

/////////////////////////////////////////////// t ////////////////// i 



o LIBSEDIV EXTENDED PRECISION DIVIDE 

///////////////////////////////////////////////»/////////////////// 
// Brief description: divide a 64 -bit quadword by an integer giving an integer quotent 
// and remainder. 

// Expected: a = integer divisor 

// b address of 64 -bit quadword 

// c = address of integer quotent 

// d = address of integer remainder 

// Result: status from lib$ediv 

/////////////////////////////////////////////////////////////////// 
Ul c_routine int routine lib$ediv (void *a, void *b, void *c, void *d) 

{ Mint *) c = (int) (*(int64 *) b / * (int *) a); 
p Mint *} d = (int) (*(int64 *) b % Mint *) a); 

jf 55 ! return 1 ; } 

fear /////////////////////////////////////////////// //////////////////// 

fi | // EMUL EXTENDED PRECISION MULTIPLY ~~ 

,fs /////////////////////////////////////////////////////////////////// 

Ss4 f J Brief description: multiply two integers giving a quadword result. 
%£3 // Expected: 

111 // a integer multiplier 

// b = integer multiplicand 

// c - integer addend 

// q = address of quadword result 

// Result: status from lib$emul 

///////////////////////////////////////////////,/////////////////// 
int routine lib$emul (int *_a, int *_b, int *c, void *_d) 
{ int 64 a = *_a; 
int 64 b = *_b? 
Hi int64 d = a * b; 

j^; if (c) 

M d += *C; 

f|J *(int64 *) _d = d; 

return 1 ; } 

/////////////////////////////////////////////// , i_j_j_j_j_j_j_j_j_j_j_j_j_j_j_j_j_j_j m 
,i EDIV EXTENDED PRECISION DIVIDE 

// / //////////////////////////////////////////// j ////////////////// / 
II Brief description: divide a 64 -bit quadword by an integer giving an integer quotent 
// and remainder. 

// Expected: a = integer divisor 

// b address of 64 -bit quadword 

// c = address of integer quotent 

// d = address of integer remainder 

// Result: status from lib$ediv 

/////////////////////////////////////////////////////////////////// 
c_routine int routine ediv (int a, int *b, int *c, int *d) 
{ return lib$ediv (&a 7 b, c, d) ; } 

/////////////////////////////////////////////////////////////////// 

n EMUL EXTENDED PRECISION MULTIP 

/////////////////////////////////////////////////////////////////// 
// Brief description: multiply two integers giving a quadword result. 
// Expected: 

// a = integer multiplier 

// b = integer multiplicand 

// c = integer addend 

// q = address of quadword result 

// Result: status from lib$emul 

/////////////////////////////////////////////////////////////////// 
coroutine int routine emul (int a, int b, int c, int *q) 



Appendix, Page 6 of 43 



1 



QIPLG Attorney Docket No:500.0028-20us 



{ return lib$emul (&a, &b, &c, q);} 

/ / / / / / / / / / / i t / / / / / // // // // // // / / / // / / / // // / / / / / i / t / / // / / // // // / / // / 



SUBX EXTENDED PRECISION SUBTRACT 

/ f /// f ////////// i ///// f //////////// f /////////// i ////////////////// / 

// Brief description: subtract two 64 -bit quadword giving quadword result 
// Expected: 

// a address of 64 -bit quadword 

// b = address of 64 -bit quadword 

// c address of 64 -bit quadword 

// Result: status from lib$subx 

///////////////////////////////////////////////'/////////////////// 
c_routine int routine lib$subx {void *a, void *b, void *c) 
{ *(int64 *) c = *(int64 *) a - *{int64 *) b; 
return 1 ; } 

///////////////////////////////////////////////>/////////////////// 



LIB$ADDX EXTENDED PRECISION ADD 



///////////////////////////////////////////////'/////////////////// 
1 1 Brief description: add two 64 -bit quadword giving quadword result 
// Expected: 

// a = address of 64 -bit quadword 

// b = address of 64 -bit quadword 

// c - address of 64 -bit quadword 

// Result: status from lib$addx 

//.///////////////////////////////// / / ///////////////////////////// / 
coroutine int routine lib$addx {void *a, void *b, void *c, void*) 
{ *{int64 *) c = *{int64 *) a + * (int64 *) b; 
return 1 ; } 

RANDOM NUMBERS 

/////////////////////////////////////////////// 
H" | \ MTH$RANDOM GET A RANDOM VALUE 

Q ////// / f //////////////// f fiii t /////////////////,////////////////// / 

'J? // Brief description: Get a random value modulus 'maxvalue 1 or if no 'maxvalue' 

12 // return a value between 0 and 1. 

5p // Expected: maxvalue = address of an integer 

'«« // Result: floating random number 

///////////////////////////////////////////////</////////////////// 

real routine mth$random {int *maxvalue) 
^ { int r = ran_random {) ; /* get random integer*/ 

jfy real one = 1; 

ff| real f ; 

if {*maxvalue) 

¥ { r %= *maxvalue; /* remainder after divide by maxvalue*/ 

Q f = (real) {r + 1); /* relative to one*/ } 

else 

{ f = (real) r; 

while (f >= one) /* make the random number less than one*/ 

Rj f /= 10 <* } 

5s7| return f; /* return floa/real*/} 

$»» # include <stdio.h> 

jfy /* randomx: 



An improved random number generation package. 

/////////////////////////////////////////////////////////////////// 
/* In addition to the standard rand{ } /srand( ) like interface, this package also has a special state info 
interface. The initstateO routine is called with a seed, an array of bytes, and a count of how many 
bytes are being passed in; this array is then initialized to contain infcrmation for random number 
generation with that much state information. Good sizes for the amount of state information are 32, 64, 
128, and 256 bytes. The state can be switched by calling the setstateO routine with the same array as 
was initiallized with initstateO . By default, the package runs with 128 bytes of state information and 
generates far better random numbers than a linear congruential generator. If the amount of state 
information is less than 32 bytes, a simple linear congruential R.N. G. is used. 

Internally, the state information is treated as an array of longs; the zeroeth element of the array is the 
type of R.N.G. being used {small integer) ; the remainder of the array is the state information for the 
R.N. G. Thus, 32 bytes of state information will give 7 longs worth of state information, which will allow 
a degree seven polynomial. (Note: the zeroeth word of state information also has some other information 
stored in it -- see setstateO for details). 

The random number generation technique is a linear feedback shift register approach, employing trinomials 
(since there are fewer terms to sum up that way) . In this approach, the least significant bit of all the 

numbers in the state table will act as a linear feedback shift register, and will have period 2^deg - 1 
{where deg is the degree of the polynomial being used, assuming that the polynomial is irreducible and 

primitive) . The higher order bits will have longer periods, since their values are also influenced by 

pseudo-random carries out of the lower bits. The total period of the generator is approximately 
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deg*(2**deg - 1); thus doubling the amount of state information has a vast influence on the period of the 
generator . 

Note: the deg*(2**deg - 1) is an approximation only good for large deg, when the period of the shift 
register is the dominant factor. With deg equal to seven, the period is actually much longer than the 
7* (2**7 - i) predicted by this formula. 

For each of the currently supported random number generators, we have a break value on the amount of state 
information (you need at least this many bytes of state info to support this random number generator) , a 
degree for the polynomial (actually a trinomial) that the R.N.G. is based on, and the separation between 
the two lower order coefficients of the trinomial.*/ 



/////// 


//////// 


/ / / 


/ / / 


/////////////////// 


#def ine 


TYPE_0 


0 


/* 


linear congruent ial*/ 


#def ine 


BREAK_0 


8 






#def ine 


DEG_0 


0 






#def ine 


SEP_0 


0 






#def ine 


TYPE 1 


1 


/* 


x**7 + x**3 + 1*/ 


#def ine 


BREAK_1 


32 






#def ine 


DEG_1 


7 






#def ine 


SEP_1 


3 






#def ine 


TYPE_2 


2 


/* 


x**15 + x + 1*/ 


#def ine 


BREAK_2 


64 






#def ine 


DEG_2 


15 






#def ine 


SEP_2 


1 






#def ine 


TYPE_3 


3 


/• 


x**31 + x**3 + 1*/ 


#def ine 


BREAK_3 


128 






#def ine 


DEG_3 


31 






#define 


SEP_3 


3 






#def ine 


TYPE_4 


4 


/* 


x**63 + x + 1*/ 


#define 


BREAK_4 


256 






#def ine 


DEG_4 


63 






#def ine 


SEP_4 


1 






/* * Array versions of 


the above 


information to make code run 


* on fact 


that TYPE_i 




.*/ 




#def ine 


MAX_TYPES 


5 


/* 


max number of types above*/ 



relies 



static int degrees [MAX_TYPES] = { DEG_0, DEG__1, DEG_2,DEG_3, DEG__4 } ; 
static int seps [MAX_TYPES] = { SEP_0, SEP_1, SEP_2,SEP_3, SEP_4}; 
/* * Initially, everything is set up as if from : 

* initstate{ 1, &randtbl, 128 ); 

* Note that this initialization takes advantage of the fact that ran_s random ( ) 

* advances the front and rear pointers 10*rand_deg times, and hence the 

* rear pointer which starts at 0 will also end up at zero; thus the zerceth 

* element of the state information, which contains info about the current 

* position of the rear pointer is just 

* MAX_TYPES* {rptr - state) + TYPE_3 == TYPE_3.*/ 
static long randtbl [DEG_3 + 1] = { TYPE_3, 

0x9a319039, 0x32d9c024, 0x9b663182, 0x5dalf342, 
0xde3b81e0, 0xdf0a6fb5, 0xfl03bc02, 0x48f340fb, 
0x7449e56b, OxbebldbbO, 0xab5c5918, 0x946554fd, 
0x8c2e680f, 0xeb3d799f, 0xbllee0b7, 0x2d436b86, 
0xda672e2a, 0xl588ca88, 0xe369735d, 0x904f35f7, 
0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc, 
0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b, 
0xf5ad9d0e, 0x8999220b, 0x27fb47b9}; 
/* * fptr and rptr are two pointers into the state info, a front and a rear 

* pointer. These two pointers are always rand_sep places aparts, as they cycle 

* cyclically through the state information. (Yes, this does mean we could get 

* away with just one pointer, but the code for random () is more efficient this 

* way) . The pointers are left positioned as they would be from the call 

* initstate( 1, randtbl, 128 ) 

* (The position of the rear pointer, rptr, is really 0 (as explained above 

* in the initialization of randtbl) because the state table pointer is set 

* to point to randtbl [1] (as explained below).*/ 
static long *fptr = fcrandtbl [SEP_3 + 1] ; 

static long *rptr = Srandtbl [1] ; 

/* * The following things are the pointer to the state information table, 

* the type of the current generator, the degree of the current polynomial 

* being used, and the separation between the two pointers. 

* Note that for efficiency of randomO, we remember the first location of 

* the state information, not the zeroeth. Hence it is valid to access 

* state [-1] , which is used to store the type of the R.N.G. 

* Also, we remember the last location, since this is more efficient than 
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* indexing every time to find the address of the last element to see if 

* the front and rear pointers have wrapped.*/ 
static long *state = &randtbl[l]; 

static int rand_type = TYPE_3; 

static int rand_deg - DEG_3; 

static int rand_sep = SEP_3; 

static long *end_ptr = &randtbl [DEG_3 + 1] ; 

/* * ran_srandom: 

* Initialize the random number generator based on the given seed. If the 

* tv P e is tne trivial no-state-information type, just remember the seed. 

* Otherwise, initializes state [] based on the given "seed" via a linear 

* congruential generator. Then, the pointers are set to known locations 

* that are exactly rand_sep places apart. Lastly, it cycles the state 

* information a given number of times to get rid of any initial dependencies 

* introduced by the L.C.R.N.G. 

* Note that the initialization of randtbl [] for default usage relies on 

* values produced by this routine.*/ 
void ran_s random (unsigned x) 

{ register int i, j; 

if (rand_type == TYPE_0) 

{ state [0] - x ; } 

else 

{ j = i; 

state [0] = x; 

for (i - 1; i < rand_deg; i++) 

{ state [i] = 1103515245 * state [i - 1] + 12345; } 

fptr = Estate [rand_sep] ; 
rptr = fcstate [0] ; 

for (i = 0; i < 10 * rand_deg; i++) 
ran_random () ; 

}} 

/* * initstate: 

* Initialize the state information in the given array of n bytes for 

* future random number generation. Based on the number of bytes we 

* are given, and the break values for the different R.N.G.'s, we choose 

* the best (largest) one we can and set things up for it. ranjs random ( ) is 

* then called to initialize the state information. 

* Note that on return from ranjsrandomO , we set state [-1] to be the type 

* multiplexed with the current value of the rear pointer; this is so 

* successive calls to initstate () won't lose this information and will 

* be able to restart with setstateO . 

* Note: the first thing we do is save the current state, if any, just like 

* setstateO so that it doesn't matter when initstate is called. 

* Returns a pointer to the old state.*/ 

char *ran_initstate (unsigned seed, /* seed for R. N. G.*/ 

char *arg_state, /* pointer to state array*/ 
int n) /* # bytes of state info*/ 

{ register char *ostate = (char *) (&state [-1] ) ; 
if (rand_type == TYPE_0 ) 

state [-1] = rand_type; 
else 

state [-1] = MAXJTYPES * (rptr - state) + rand_type; 
if (n < BREAK_1) 

{ if (n < BREAK_0) 

{ fprintf (stderr, 

"initstate: not enough state (%d bytes) with which to do jack; ignored." 
n); 

return 0; } 
rand_type = TYPE_0; 
rand_deg ~ DEG_0; 
rand_sep = SEP_0; } 

else 

{ if (n < BREAK_2) 

{ rand_type = TYPE_JL; 

rand_deg = DEG_1; 
rand_sep = SEP_1; } 

else 

{ if (n < BREAK_3) 

{ rand_type = TYPE_2; 

rand_deg = DEG_2 ; 
rand_sep = SEP_2 ; 
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else 

{ if (n < BREAK_4) 

{ rand_type = TYPE_3; 

rand_deg = DEG__3 ; 
rand_sep = SEP_3 ; 

} 

else 

{ rand_type = TYPE_4; 

rand_deg = DEG_4; 
rand_sep = SEP_4 ; 

}}}} 

state = &{((long *) arg_state) [1] ) ; /* first location*/ 

end_j>tr = Estate [rand_deg] ; /* must set end_ptr before ran_s random*/ 

ran_s random (seed) ; 

if (rand_type == TYPE_0) 

state [-1] = rand_type; 
else 

state [-1] = MAX_TYPES * (rptr - state) + rand_type; 
return {ostate) ; } 
/* * setstate: 

* Restore the state from the given state array. 

* Note: it is important that we also remember the locations of the pointers 

* in the current state information, and restore the locations of the pointers 

* from the old state information. This is done by multiplexing the pointer 

* location into the zeroeth word of the state information. 

* Note that due to the order in which things are done, it is OK to call 
| - * setstate 0 with the same state as the current state. 

4 IT * Returns a pointer to the old state information.*/ 
0 char *ran_setstate (char *arg_state) 
q { register long *new_state = (long *) arg_state; 
?S register int type = new_state[0] % MAX_TYPES; 

1M register int rear = new_state [0] / MAXJTYPES; 

*fj char *ostate = (char *) (&state [-1] ) ; 

!( f| if (rand_type == TYPE_0) 

** state [-1] = rand_type; 

else 

||1 state [-1] = MAX_TYPES * (rptr - state) + rand_type ; 

switch (type) 
* { case TYPE_0: 

Q case TYPE 1 : 



y 



case TYPE_2: 
case TYPE_3: 
f& case TYPE_4 : 

|lj rand_type = type; 

^ rand_deg = degrees [type] ; 

rand_sep = seps [type] ; 
fU break; 

default : 

f print f (stderr, 

"setstate: state info has been munged; not changed."); } 
state = &new_state [1] ; 
if (rand_type != TYPEJ)) 

{ rptr - &state [rear] ; 

fptr » &state[ (rear + rand_sep) % rand_deg] ; } 
end_ptr = &state [rand_deg] ; /* set endjptr too*/ 
return (ostate) ; } 
/* * random: 

* If we are using the trivial TYPE_0 R.N.G., just do the old linear 

* congruential bit. Otherwise, we do our fancy trinomial stuff, which is the 

* same in all ther other cases due to all the global variables that have been 

* set up. The basic operation is to add the number at the rear pointer into 

* the one at the front pointer. Then both pointers are advanced to the next 

* location cyclically in the table. The value returned is the sum generated, 

* reduced to 31 bits by throwing away the "least random" low bit. 

* Note: the code takes advantage of the fact that both the front and 

* rear pointers can't wrap on the same call by not testing the rear 

* pointer if the front one has wrapped. 

* Returns a 31-bit random number.*/ 
long routine ranjrandom () 

{ long i; 
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if {rand_type == TYPE_0) 



{ 

else 
{ 



State [0] = (State [0] * 1103515245 + 12345) & 0x7fffffff; 



m 

S 



*fptr += *rptr; 

i = <*fptr » 1) & 0x7fffffff ; /* chucking least random bit*/ 
if (++fptr >= end_ptr) 
{ fptr = state; 

++rptr; } 

else 

{ if (++rptr >= end_ptr) 

rptr - state; 

} } 

return ( i ) ; } 
static unsigned char const ascii_to_ebcdic [] = 
{ 0, 01, 02, 03, 067, 055, 056, 057, 

026, 05, 045, 013, 014, 015, 016, 017, 

020, 021, 022, 023, 074, 075, 062, 046, 



030, 031, 077, 047, 034, 035 

0100, 0117, 0177, 0173, 0133 
0115, 0135, 0134, 0116, 0153 
0360, 0361, 0362, 0363, 0364 
0370, 0371, 0172, 0136, 0114 
0174, 0301, 0302, 0303, 0304 
0310, 0311, 0321, 0322, 0323 
0327, 0330, 0331, 0342, 0343 
0347, 0350, 0351, 0112, 0340 
0171, 0201, 0202, 0203, 0204 
0210, 0211, 0221, 0222, 0223 
0227, 0230, 0231, 0242, 0243 
0247, 0250, 0251, 0300, 0152 
040, 041, 042, 043, 044, 025 
050, 051, 052, 053, 054, Oil 
060, 061, 032, 063, 064, 065 
070, 071, 072, 073, 04, 024, 

0101, 0102, 0103, 0104, 0105 
0111, 0121, 0122, 0123, 0124 
0130, 0131, 0142, 0143, 0144 
0150, 0151, 0160, 0161, 0162 
0166, 0167, 0170, 0200, 0212 
0216, 0217, 0220, 0232, 0233 
0237, 0240, 0252, 0253, 0254 
0260, 0261, 0262, 0263, 0264 
0270, 0271, 0272, 0273, 0274 
0312, 0313, 0314, 0315, 0316 
0334, 0335, 0336, 0337, 0352 



036, 037, 

0154, 0120, 0175, 

0140, 0113, 0141, 

0365, 0366, 0367, 

0176, 0156, 0157, 

0305, 0306, 0307, 

0324, 0325, 0326, 

0344, 0345, 0346, 

0132, 0137, 0155, 

0205, 0206, 0207, 

0224, 0225, 0226, 

0244, 0245, 0246, 

0320, 0241, 07, 
06, 027, 
012, 033, 
066, 010, 
076, 0341, 

0106, 0107, 0110, 

0125, 0126, 0127, 

0145, 0146, 0147, 

0163, 0164, 0165, 

0213, 0214, 0215, 

0234, 0235, 0236, 

0255, 0256, 0257, 

0265, 0266, 0267, 

0275, 0276, 0277, 

0317, 0332, 0333, 

0353, 0354, 0355, 

0375, 0376, 0377}; 



0356, 0357, 0372, 0373, 0374 
static unsigned char const ebcdic_to_ascii [] 
{ 0, 01, 02, 03, 0234, 011, 0206, 0177, 

0227, 0215, 0216, 013, 014, 015, 016, 017, 
020, 021, 022, 023, 0235, 0205, 010, 0207, 
030, 031, 0222, 0217, 034, 035, 036, 037, 
0200, 0201, 0202, 0203, 0204, 012, 027, 033, 
0210, 0211, 0212, 0213, 0214, 05, 06, 07, 
0220, 0221, 026, 0223, 0224, 0225, 0226, 04, 
0230, 0231, 0232, 0233, 024, 025, 0236, 032, 
040, 0240, 0241, 0242, 0243, 0244, 0245, 0246, 
0247, 0250, 0133, 056, 074, 050, 053, 041, 
046, 0251, 0252, 0253, 0254, 0255, 0256, 0257, 
0260, 0261, 0135, 044, 052, 051, 073, 0136, 
055, 057, 0262, 0263, 0264, 0265, 0266, 0267, 
0270, 0271, 0174, 054, 045, 0137, 076, 077, 
0272, 0273, 0274, 0275, 0276, 0277, 0300, 0301, 

0302, 0140, 072, 043, 0100, 047, 075, 042, 

0303, 0141, 0142, 0143, 0144, 0145, 0146, 0147, 
0150, 0151, 0304, 0305, 0306, 0307, 0310, 0311, 
0312, 0152, 0153, 0154, 0155, 0156, 0157, 0160, 
0161, 0162, 0313, 0314, 0315, 0316, 0317, 0320, 
0321, 0176, 0163, 0164, 0165, 0166, 0167, 0170, 
0171, 0172, 0322, 0323, 0324, 0325, 0326, 0327, 
0330, 0331, 0332, 0333, 0334, 0335, 0336, 0337, 
0340, 0341, 0342, 0343, 0344, 0345, 0346, 0347, 
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0173, 0101, 0102, 0103, 0104, 0105, 0106, 0107, 
0110, 0111, 0350, 0351, 0352, 0353, 0354, 0355, 
0175, 0112, 0113, 0114, 0115, 0116, 0117, 0120, 
0121, 0122, 0356, 0357, 0360, 0361, 0362, 0363, 
0134, 0237, 0123, 0124, 0125, 0126, 0127, 0130, 
0131, 0132, 0364, 0365, 0366, 0367, 0370, 0371, 
060, 061, 062, 063, 064, 065, 066, 067, 
070, 071, 0372, 0373, 0374, 0375, 0376, 0377}; 
/* * lib$tra_asc_ebc translate ASCII to EBCDIC 
**/ 

int routine lib$tra_asc_ebc (DESC_S *source, DESC_S *dest) 
{ unsigned char *asc, *ebc; 
int slen, dlen; 
int i; 

split$desc (source, &slen, &asc) ; 
split$desc (dest, fcdlen, &ebc) ; 
for (i = 0; i < slen && i < dlen; i++) 
ebc[i] = ascii_to_ebcdic [asc [i] ] ; 

return SS$_NORMAL; } 
/* * lib$tra_ebc_asc translate EBCDIC to ASCII 
**/ 

int routine lib$tra_ebc_asc (DESC_S *source, DESC_S *dest) 
{ unsigned char *asc, *ebc; 
int slen f dlen; 
int i; 

split$desc (source, fcslen, &ebc) ; 
split$desc (dest, &dlen, &asc) ; 
for (i = 0; i < slen && i < dlen; i++) 
CI asc [i] = ebcdic_to_ascii [ebc [i] ] ; 

return SS$_NORMAL;} 
/* * Convert packed decimal to leading separate 



yp * string: pointer to output ascii string buffer 

if* * packed: pointer to input packed decimal string*/ 

Jft int routine cvtps(char *string, char *packed) 

{ int leading = 1; 
iQ int negative = 0; 

char buf [32+1] ; 
char *s = buf; 
int i , c ; 

for (i - 0; i < 32; i++) { if (i&D /* odd is lower*/ 
c = packed [i/2] & OxOf; 
else /* even*/ 
fy c = (packed [i/2] »4) & OxOf; 

iws, if (leading && c == 0) 

M continue; /* skip leading 0's*/ 

leading =0; /* nolonger leading*/ 
if (c == 10 | | c == 12 | | c == 14 | | c == 15) 

break; /* plus*/ 

else if (c == 11 | | c == 13) { negative = 1; 

break; } 
else if (c > 9) 

break; 
*s++ = (char) (c + ■ 0 ■ ) ; 
} 

*S++ = ' \0 ' ; 
if (negative) 

*string++ ■ * -■; /* no unary '+'*/ 
for (i = 0; buf [i] ; i++) 

*string++ = buf [i] ; 
if (i == 0) 

*string++ = '0'; /* at least one zero*/ 
*string++ = ! \0' ; 
return SUCCESS ; } 
/* * Convert leading separate to packed decimal 
**/ 

int routine cvtsp(char *packed, char *string, int len) 
{ int negative = 0; 

char buf [32+1] ; 

char *s = buf; 
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int i; 

for (i = 0; i < len && s < buf+32; { if (string [i] <= 1 ' ) 

continue; 
else if (string [i] ' .') 

cont inue ; 
else if (string [i] == 1 -') 

negative = 1; 
else if (string [i] == '+') 

negative = 0; 
else if (string [ij >= '0' && string [i] <= '9') 

*s++ = string [i] ; 

} 

*S++ = ' \0« ; 

len = strlen(buf); 

/* * if first zero padding required*/ 

if ((len&l) == 0) { for (i = 0; buf [i] ; i++) ; 

for {; i >= 0; --i) 
buf [i+1] = buf [i] ; 

buf[0] = ' 0'; 

len = strlen(buf ) ; 

for (i - 0; i < len; i++) { if (i&D /* odd is lower*/ 

packed [i/2] |= (buf [i] &0x0f ) ; 
else 

packed[i/23 - (char) ( (buf [i] &0x0f) «4) ; 

} 

packed [i/2] |= negative ? 13 : 12; 
^ return (i/2)+l; /* return length of packed decimal string*/} 

LJ /* * convert packed decimal to real 
O * 

« * packed: pointer to input packed decimal string 
*5f * result: pointer to output result real*/ 
SB int routine cvtpdr(real *result, char *packed, int digits) 
;,f| { char string [32+1] ; /* maximum packed decimal length*/ 
ga s int i, n; 

/* * Convert packed decimal to leading separate 

*/ 

cvtps (string, packed); 

/* * se t packed decimal precision to 'digits' fractional digits*/ 
for (i = 0; string[i]; i++) ; 
for (n = 0; n < digits && i >= 0; --i, n++) 

string [i+1] = string [i] ; 
string [i+1] = ' . 1 ; 

/* * Convert ascii string to real*/ 
return ator ( result, string);} 
%j?r j ± + c OIlver t real to packed decimal 



m 



ri 



* packed: pointer to output packed decimal string 

* digits: sizeof output buffer (as packed decimal digits) 

* val: pointer to real value to convert*/ 

int routine cvtrpd (char *packed, int digits, int precision, real *val) 
{ real info info; 

char buf [32+1] ; /* maximum size of a packed decimal field*/ 

char pbuf [32+1] ; /* packed buffer*/ 

int p_digits; 

DESCJ3 desc; 

int i; 

set$desc (&desc, sizeof (buf), buf ) ; 

/* * Convert real to ascii string*/ 

rtoa(val, &desc, precision, &info f 0) ; 

/* * Convert ascii string to packed decimal*/ 

p_digits = cvtsp(pbuf, buf, 0) ; 

digits = (digits/2) +1; /* convert digits to bytes + sign*/ 

/* * fill packed decimal string with required 

* leading packed zeros.*/ 
for (i = p_digits; i < digits; i++) 

*packed++ = ' \0 ■ ; 
/* * now fill with converted packed number*/ 
for (i = 0; i < p_digits; i++) 

*packed++ = pbuf[i]; 



Appendix, Page 13 of 43 



QIPLG Attorney Docket No:500.0028-20us 



return SUCCESS;} 
/* * cvtlp -- convert long to packed decimal 
**/ 

int routine cvtlp (int *expr, int digits, char *buf) 
{ real r = *expr; 

return cvtrpd(buf / digits, 0, &r) ; } 
/* * cvtpl -- convert packed decimal to long 
* */ 

int routine cvtpl (int digits, char *buf, int *result) 
{ int stat; 
real r; 

stat = cvtpdr(&r, buf, digits); 
*result = (int)r.ip; 
return stat ; } 

DEBUG DUMP HEX BYTES ,,,,,,, 

///////////////////////////////////// /// ^ / 11 1 ' ' 1 ' 1 ' ' 1 1 f 1 / f 1 1 f 1 ' 1 1 

„ PHEXB DU MP HEX BYTES OF A MEMORY LOCATION ,,,,,, ' 

//////////// 7 ////////////////////// //// / / / ft * > i • t i i i 'i ' i ' i ' ' ' > 111 ' ' 

// Brief description: This function is a debug function. It will dump in hexadecimal 

// format with ascii characters at the right, a block a memory 

// starting at address 'a' for bytes 'n' . Unprintable character 

// hex values are displayed but their ascii value is 

// This function is called when debugging the math routines 

// at a very low level. It is called to dump the contents of 

// 'real' 64-bit integers (and their 128-bit intermediate values) 

// to check that the hex values are indeed correct. 

// Any block of memory may be dumped with this function; not 

// just 'real' 64 -bit (or intermediate 12 8 -bit) values. I have 

// used this function to dump »pcode' blocks of memory when I 

// was not sure what was being generated. 

// Expected: a = the start address (any type) 

// n = number of bytes to dump. 

// / Result: / return , 

void routine dhexb (void *a, int n) 

{ unsigned char *s - (unsigned char *) a; 

char msg [256] ; 

char buf [32] ; 

int sn = n; 

int i; 

for {; n > 0; n -= i) 

{ sprintf (msg, 11 %04d » , sn - n) ; 

for (i = 0; i < 16 && i < n; i++) 
{ sprintf (buf, "%02X " , s [ij ) ; 

strcat (msg, buf) ; } 
for (; i < 16; i++) 

strcat (msg, 11 ") ; 
strcat (msg, " "); 

for (i = 0; i < 16 && i < n; i++, s++) 
{ sprintf (buf, »%c», *s< 1 1 j | *s> 1 - 1 ? '.' : * s I ; 

strcat (msg, buf) ; } 
logf ("%s\n" , msg) ; 

, , , )) / / y /////////////////////////// / / ////// / / / ' ' * 1 1 1 ' 1 ' ' 1 1 ' 1 11 ' 1 ' ' ' 

// Brief description: add 2 reals and place result in a third 
// Expected: 1st arg = unused 

// a = the address of a quadword (integer pointer) 

// b = the address of a quadword (integer pointer) 

// c = the address of a quadword (integer pointer) 

// 

Z 7 / / / R ? S / U Vy / / C / / / / / / / = / a / / ?'/ /////////// / / / / / ///////// / / / / / / / / / / " " ' ' ' ' ' ' 
void routine addm (int, int *a, int *b, int *c) 
{ *(int64 *) c = *(int64 *) a + *(int64 *) b;} 



MORE MATH AND COMPARE FUNCTIONS 

f f l l 1 I ! t t I ! I I I I l l t t t t I t f I I I * I I ' ' > 1 1 1 

SUBM SUBTRACT 2 QUADWORDS AND PLACE RESULT IN A THIRD 



, / 7 j > , / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' ' i ' ' ' * 1 ' ! < > ' 1 1 1 ' 1 ' 1 1 ' ' 1 1 1 1 ' ' ' 



", " /'"/ /'/ /" '/ '/ /' / 1 i f T '/ 7-1 n ' / ( / / / / / 7 i / / / / / / / / / / / / / / / / / ////////////////'''> 

// Brief description: subtract 2 reals and place result in a third 
// Expected: 1st arg = unused 

II a = the address of a quadword (integer pointer) 
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// b = the address of a quadword (integer pointer) 

// c = the address of a quadword (integer pointer) 

// 

// Result: c = b - a; 

///////////////////////////////////////////////'////// / ' ' f / { * ' / 1 1 * 1 
void routine subm (int, int *a, int *b, int *c) 
{ *(int64 *) c = *{int64 *) b - * (int64 *) a;} 

/ ////////////////////////////////////// / til / / / / t ttt//l/t//It/ftttff 

. IS NEGATIVE IS A REAL NEG ATIVE _ 

////////////////////////////// / / ////////////// / / // // // // // // // // // / 
// Brief description: test the sign bit of a real 
// Expected: a = pointer to a real 
// Result: true if negative 

J////////////////////////////////////////////////////' ft 'i 1 1 11 1 1 1 1 1 
inline int routine i s_nega tive_real (real * a) 
{ if (a->ip) 

return a->ip < 0; 
else 

return a->fp < 0;} 

/////////////// / t / // // // / / //////////// / t / i l f l t f t f f f ' f f f t / / / / / / / / / / / 

u IS_NEGATIVE IS A 64-BIT INTEGER ARR AY NEGATIVE 

////,/////////////////// /////////////////////// //////////////////// 
// Brief description: test the sign bit of a 64-bit integer array 
// Expected: a = pointer to start of a 64-bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: true if negative 

//////////////////////////////' t ' t / / ////////////// / / / / / / // // // // // ' 
inline int routine is_negative (void *a, int n) 

{ return (int) ( ( ( (uint64 *> a) [n - 1] » 63) & 1); /* shift the 64-bit integer down*/} ,,,,,,,, 
/ //////////////////////////////////// f f ' f f t f f f t i t / i / i i / / / / / / / / / L ' 1 1 



IS A REAL LESS THAN ANOTHER 



///////////////////////////////////// ////////// //////// / / / / / / ' ' / / ' 1 
II Brief description: is a < b 
// Expected: a = pointer to a real 
// b ■ pointer to a real 

/'/ / z**? 8 / 11 /^" / / / r f 8U /" fc / f /^/H ///////////////////// / / / / ' ' ' / // // // // // // / / / / / / t 
inline int lss (real * a, real * b) 
{ if (a->ip < b->ip) 

return 1; 
if (a->ip > b->ip) 

return 0; 

/ / I^juLlj. f</ J ^Pf ]f /////////'/ / / i f L t t f f f t / //////////// f f f f i f f i t ' i 1 i 1 1 
u NEGATE A 64-BIT INTEG ER ARRAY . 

//////////J//////////////////////// ''' f f f ' f f t f f /////// / / / / ' ' ' ' ' ' ' ' ' 

// Brief description: negate a 64-bit integer array (make negative positive or 
// positive negative) 

// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

/ / R ? S / U / t / / / / V ?^/ / // // // // // // // // // // // // // // // // / / / / / / / ' / / // // // // / / 
inline void routine negate (void *a, int n) 
{ int i; 

({int64 *) a) [03 = -{(int64 *) a) [0] ; /* two's compliment*/ 
for (i = 1; i < n; i++) 

((int64 *) a) [i] = ~{{int64 *) a) [i] ; /* one's complement*/} ,,,,,,,,,,,,, 
/ / / / / /,///// / //////////////// / /////// t ' t fit i t / / / i i / t i / / / i / / / / / / / 

COPY A 64-BIT INTEGER ARRAY 



////,/// * ///////////////////// ////////// t' t / ////////////////'/// 
// Brief description: inline copy a number of 64 -bit integers from 'b' to 'a' 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

// n length of array (number of 64-bit integers) 

/ ^ / /w^/V / / / V ?V*/ / // // // // // // // // // / / / // // // / / ///////ft/ J / / // // // // / / 
inline void routine cpy (void *a, void *b, int n) 
{ int i; 

for (i - 0; i < n; i++) 

((int64 *) a) [i] = ( (int 6 4 *) b)[i];} 

/ l l I ! I f I I I J I / / / / / / I / / / / / / / / / / / / / / / / / / / / / / f f t I t t t t t / / / ////// // / / / / / / 



/, SHIFTUP BY ONE BIT A 64-BIT INTEGER ARRAY 

////////// / /////// f / /////// 7 / / / / / / / / / / / / / / / / / / / ////////// f f f t t t f f f ' 
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// Brief description: shiftup by one bit a 64-bit integers array 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: void 

//////////////////////////// / / ft / // // / / ft it / / / / > i ttiititiitttiiitii 
inline void routine shiftup (void *a, int n) 
{ int c, carry = 0; 
int i; 

for (i = 0/ i < n; i++) 

{ c = (int) (((uint64 *) a) [i] >> 63) & 1; /* shift the 64-bit integer down*/ 

( (int64 *) a) [i] «= 1; 
((int64 *) a) [i] |= carry? 
carry = c; 

i , , )\ , I ( i i t j i i i / i t ! i t t i i t ! I I I ! i t f i f i i i f i t i i i t t i >' f i * f t * i 1 1 1 1 ' 1 1 1 f 1 f 1 
u SHIFTDOWN BY ONE BIT A 64-BIT INTEGE R ARRAY 

////////////////////////////////////////// ////// t ! i ! f t t t t l f ) t i f I f f t 

f I Brief description: shiftdown by one bit a 64 -bit integers array 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

i^t t / R ? S / U /V / ( / v / 0: / d / tiitifttttftfttttft 1 1 til / / / / / / / / / ///////// / ft/ffttfft 

inline void routine shiftdown (void *a, int n) 
{ int 64 c, carry = 0; 
int i; 

for (i = n - 1; i >= 0; --i) 

{ c = (int) (((uint64 *) a) [i] ) & 1; 

((uint64 *) a) [i] »= 1; /* unsigned shift down*/ 
((uint64 *) a) [i] |= (carry << 63); 
carry = c; 



h LSS COMPARE FUNCTIONS FOR AN ARR AY OF 64-BIT INTEGERS 

/,//,//////,//////////// //// / / / / f f /// f t f ' t f ' f f f ' i i f f i f f * i i i f i i i i 1 f * 

H LSS LESS THAN COMPARE AN ARRAY OF 64-BIT INTEGERS ,,,,,,,, y ,,//// / 

/////////////////////f///tt/tft//ft / i it f f i i i i i i tfftffffffttttftftf/ 

II Expected: 

// a address of 64 -bit integer array 

// b address of 64-bit integer array 

// n number of 64 -bit integers 

// Result: TRUE if (a < b) else FALSE; 

/ / / / / / / ////// / / // /t//f//t/t////ff//ffy f ft/ t i / / /f/t/ftstf/ft//t///ft 
inline int routine lss (void *a, void *b, int n) 
{ int i; 

if (((int64 *) a) [n - 1] 1= ( (int64 *) b) [n - 13) 

return {(int64 *) a) [n - 1] < ( (int64 *) b) [n - 1] ; /* signed conpare*/ 
for (i = n - 1; i >= 0; --i) 

{ if (({int64 *) a) [i] I - { (int 64 *) b) [i] ) 

return ( (uint64 *) a) [i] < ( (uint64 *) b) [i] ; /* unsigned compare*/ } 
return 0; /* false equal*/} 

f f l I f I i / I t l t I f f I f i I f I l t t I f f I Iff I l f ! / / / / / // // // /f/fftfftfftft/ff/fff 

„ LEG LESS THAN OR EQU AL COMPARE TWO REALS 

y /// y ////////////////////////////////////////// / /////////////////// 

// Expected: 

// a address of real 

// b address of real 

// Result: TRUE if (a <= b) else FALSE ; ,,,,,,,,,,,,,,,,,,,, 
f////////////////////t//tff//fttf///ftttft/tft/f/tf/ffttt//t//ff/ff 

inline int routine leq (real * a, real * b) 

{ if (a->ip <= b->ip) 

return 1; 

if (a->ip > b->ip) 

return 0; 

/ / J^ff A^F ^7 / b 7 ////////////////// / > ' >>><>>! > > i i I t i t if / / / / / / / / / / L 

o LEQ LESS THAN OR EQUAL COMPARE AN A RRAY OF 64-BIT INTEGERS 

//////////////////// / ////////////// //// n f / / / / / /////////////////''' 
// Expected: 

// a address of 64 -bit integer array 

// b address of 64 -bit integer array 

// n number of 64 -bit integers 

// Result: TRUE if (a <= b) else FALSE; , t , i i i i t t 

'/ ,/,/////,////////// t / t t / t / t // f /// t t /// t / t f /// t // t // t / t f ' f t ' t f '''' ' 
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inline int routine leq (void *a, void *b, int n) 
{ int i ; 

if {({int64 *) a) [n - 1] 1= ( (int64 *) b) [n - 1]) 

return ( <int64 *) a) [n - 1] < ((int64 *) b) [n - 1]; /* signed compare*/ 
for (i = n - 1; i >= 0; --i) 

{ if {{(int64 *) a) [i] != ( (int64 *) b) [i] ) 

return { (uint64 *) a) [i] < ( {uint64 *) b) [i] ; /* unsigned compare*/ } 

return 1* /* true equal*/} 

//,,/////////////////////////////////// /////''/''/ ! tf 11 ' i ' ! 1 ' f 1 1 1 1 1 

. GEQ GREATER THAN OR EQUAL COMP ARE AN ARRAY OF 64-BIT INTEGERS 

//////////////////////////////// / ////////////// : /////////////////// 
// Expected; 

// a address of 64 -bit integer array 

// b address of 64 -bit integer array 

// n number of 64 -bit integers 

// Result: TRUE if (a >= b) else FALSE; ,,,,,,,,,,,,,, 
J/////,//////////////////// ////////////////// ////// ' / / / ' ' ' > 1 1 1 1 ' 
inline int routine geq (void *a, void *b, int n) 
{ int i; 

if (({int64 *) a) [n - 1] 1= {(int64 *) b) [n - 1]) 

return ( (int 64 *) a) [n - 1] > ( (int 64 *) b) [n - 1]; /* signed compare*/ 
for (i = n - 1; i >= 0; --i) 

{ if (((int64 *) a) [i] != ((int64 *) b) [i] ) 

return ( (uint64 *> a) [i] > ( (uint64 *) b) [ij ; /* unsigned compare*/ } 

return 1- /* true eoual*/} 

/ / 7 / 7 j _i 1 ,,,,,,, t i i , i i i i ii i i > ' i i i ' ><>'>><' >'!>>'■><! 1 > 1 > '> { 1 1 f 1 1 1 ' 

u SET TO ZERO A 64-BIT INTEGER ARRAY 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i / 7 i i i i i i i i i / ' i i i i i i i i i i i i i i i i i i 1 

fl Brief description: set to zero a 64-bit integer array 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit * integer array 

// (or a real) . 

// n length of array (number of 64-bit integers) 

{ f t i / R ? s / U /V/ / / v / 01 / d / /////////////////////////////////' ft I** 1 1 1 1 1 1 1 1 1 / / / / / 

inline void routine set_zero (void *a, int n) 
{ int i; 

for (i = 0; i < n; i++) 

/ / / / ) { tjL±l LJ. [/l3 ' 1 '°< t i ' ' i i < > < ' ' ' ' * > > < * > > > ! > 1 ! 1 11 11 ( 1 1 { ' 1 ' 1 1 1 1 ' ' 1 1 ' 1 ' 
n TEST IF A 64-BIT INTEGER ARRAY IS ZERO ,,,,,,,,,, 

7 7 //////////// 1 t t t i i i i i f i i f t i i ( f i i ( i i t i * * i * ' i ' i t i t l i i i l i i ' * ' i ' ' i 1 ' J 

II Brief description: test if a 64 -bit integer array is zero 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

! , f t / fTf^i i i /7V 1 / ///////// y /////////////// ' ' /////// / ' ' ' ' / ' ' ' ' ' ' ' ' ' ' ' 1 1 
inline int routine is_zero (void *a, int n) 
{ int i; 

for (i = 0; i < n; i++) 
if {((int64 *) a) [i] ) 
return 0; 

, , fyyy, 1 // /,//,/// 1 , i<< i >,"><"<<<<<<<<<< <<<<"< << 



ADD64 ADD 2 64-BIT INTEGERS AND PLACE RESULT IN A THIRD / 

//////////////////////// / i ! I l I I I * I * I I < * * ' ' 'I ' ft I I * * * I * * * I * I > >' '> 1 > 1 

H Brief description: add 2 reals and place result in a third 
// Expected: 

// a = the address of a 64 -bit real 

// b = the address of a 64 -bit real 

// c = the address of a 64-bit real 

// 

// ^ Result: f f t f f t t f ^ t ^^ t , , , , , , , , t , , , f , g , § , , , t , f f t , t , , , , , t t t t , t t t t t t i s 
void routine add64 (real * a, real * b, real * c) 
f *(int64 *) c = *{int64 *) a + *(int64 *) b;} 

//////////////////////// / / // / / / / / titiiffiiiiiiiiii 



i i i i i t i i i i i / i i i i i 



„ SUB64 SUBTRACT 2 64-BIT INTEGERS AND PLACE RESULT IN A THIRD 



////////////////////////////////////// // 11 f ' 11 1 ' ' 1 1 ' ' ' 7 ' ' ' ' ' ' ' ' ' ' ' ' 

// Brief description: subtract 2 reals and place result in a thxrd 

// Expected: 

// a = the address of a 64-bit integer (void pointer) 

// b = the address of a 64-bit integer (void pointer) 

// c = the address of a 64-bit integer (void pointer) 
// 
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1 1 Result * C ~ a — b; 

/////////////////// / / '////////////////// / > //////// / / / / / / / / / / ' ' ' / ' ' ' ' 
void routine sub64 (void *a, void *b f void *c) 
{ *(int64 *) c = *(int64 *) a - *{int64 *) b;} 



/////////////////////////// 



i t t / i t t / i i t i i i t t i i / i i i i t t i f i i i i i t t i t i i * I 



o ADD A 64-BIT INTEGER ARRAY (LENGTH OF TWO) _ 

1 /// //////////////////////////////////////////// ' //////////////////' 
// Brief description: inline add a length of two 64 -bit integers a = a + b 
// slightly faster for a 'normal' add of a 128-bit integer 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64-bit integer array 
// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

/ / /*? S / U /V / / / V / 01 / d / //////// //////////////////////'// ' ' / / / / / / / / / / / / / / ' ' ' ' 
inline void routine add2 (void *a, void *b) 
{ uint64 was; 

int64 carry = 0; 

was - {(uint64 *) a) [0] ; 

({uint64 *) a) [0] +« ( (uint64 *) b) [0] ; 

carry = { ( (uint64 *) a) [0] < was); 

( (uint64 *) a) [1] += ( (uint64 *) b) [1] + carry;} ,,,,,,,,.,,/////,/ y //////// / 
h SUBTRAC T A 64-BIT INTEGER ARRAY (LENGTH OF TWO) _____ — 

1 1 Brief description: inline subtract a length of two 64-bit integers a = a - b 

// slightly faster for a 'normal' subtract of a 128-bit integer 

// (the pointers can point at an integer array or a real) 
s „ // Expected: a = pointer to start of a 64 -bit integer array 
\ If (or a real) . 

Q // b pointer to start of a 64-bit integer array 

G // * or a real ) * 

*J3 /( / fTY^f / / 7J^i /////////////////////////////////'/////////////////// 
..ft inline void routine sub2 (void *a, void *b) 
^ { uint64 was; 
%0 int64 carry = 0; 

ps was = ( (uint64 *) a) [0] ; 

12 ((uint64 *) a) [0] -= ( (uint64 *) b) [0] ; 

carry = ( ( (uint64 *) a) [0] > was) ; 

(<uint64 *) a) [1] -= ((uint64 *) b) [1] + carry;} ,,,,,,,,,,,,/,,,///////////// 



. SCALEUP A 1 28-BIT INTE GER BY SCALE _ 

///////////////////////////// ////////_/_/ ^ / / ' / ' 1 1 ' i t i t ' i i i i i 1 



//////// 

II Brief description: Multiply the 128-bit integer by SCALE 
H 8 // Expected: 

ni If _a = address of 64-bit integer array (_a * _b) 

O i f i / / R f ^/V / w 1 ? ////////////////////////////'/'''/' ////// / / / / / / / ' ' ' / ' ' ' 
iJV inline void routine scaleup (void *_a) 
ty { int64 a[] = { 0, 0 }; 

int64 b[] = { SCALE, 0 }; 

int64 c[3 = { 0, 0 }; 

int na = 0; 

int i; 

Cpy (a, J3l, 2) ; 

if (is_negative (a, 2)) 

{ /* negative*/ 

negate (a, 2); /* make positive*/ 
na = 1; } 
if (lss (a, b, 2)) 

{ /* swap*/ 

int 64 temp [23 ; 
cpy (temp, b, 2) ; 
cpy (b, a, 2) ; 
cpy (a, temp, 2) ; } 
for (i = 0; i < 64 && <b[0] || b[l]); i++) 
{ if (b[0] & 1) 

add2 (c, a) ; 
shiftdown (b, 2) ; 
shiftup (a, 2); } 
if (na == 1) 

{ negate (c, 2); /* make positive*/ } 

, / ? p / y / ( -?'/ C i'i V/l /,///,//////////////////////// ' ' / ' / ' ' ' ' ' ' ' ' ' ' / 1 1 1 1 1 11 1 
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SUBTRACT A 64-BIT INTEGER ARRAY 



///////////// 7 ///////////// ///////////////// ~ / < ' i * ' ' ' ' 1 ' ' ' ' ' 1 ' ' 1 1 ' 
II Brief description: inline subtract a number of 64-bit integers a = a - b 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

// n length of array (number of 64-bit integers) 

i 1 , / fVf^i i i / V / 01 / d / //////////////////////////////// ' //////// / / / / / / / / / / / / 
inline void routine sub (void *a, void *b, int n) 
{ uint64 was; 

int 6 4 carry = 0; 
int i; 

for (i = 0; i < n; i++) 

{ was = ( (uint64 *) a) [i] ; 

(<uint64 *) a) [i] ( (uint64 *) b) [i] + carry; 
carry = ( ( (uint64 *) a) [i] > was) t 

, , , )} , i , t , , i , i , t I ! I i I I l I I i t i i t i i i i i i i i t i > i * i ' > i ' < > > ' 1 1 1 1 1 1 1 ' 1 { 1 ' ( 1 ' 

» ADD A 64-BIT INTEGER ARRAY ,,,,,,, 

//,/,////////////////////////////'/// / / / / i I''' ~ ' 1 ' ^ ' 
// Brief description: inline add a number of 64-bit integers a = a + b 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

!i / fTf^'i i / / V w d / ///////////////////////' ' / ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' 

U inline void routine add (void *a, void *b, int n) 
»f| { uint64 was; 

int64 carry = 0; 
int i; 

*|J for (i = 0; i < n; i++) 

fxi { was = ( (uint64 *) a) [i] ; 

^ ((uint64 *) a) Ei3 += ( (uint64 *) b) [ij + carry; 

59 carry = { ( (uint64 *) a) [i] < was); 

* fJ-L LU ! t ! ! I i t t < f I I * > ! * ' ' ' * I ' 1 * ' 1 1 { 1 1 ' ' ' 1 ' 1 ' ' ' ' ' ' ' ' 1 ' ' ' ' ' ' ' ' ' 1 ' ' ' ' ' ' ' 

O 1 , SCALED OWN A 128-BIT INTEGER BY SCALE ,,,,,,,,,,,, 

II Expected: 

}«& // _ a = address of 64-bit integer array 

W , ///,/////////////////////////////'//""""""""" 

m inline void routine scaledown (void *_a, int n) 
^1 { uint64 a [] = { 0, 0, 0, 0 }; 
iw uint64 r[3 = { 0, 0 }; 

uint64 scale [] = { SCALE, 0, 0 f 0 }; /* SCALE must be less than an mt64*/ 

uint64 b[] = { SCALE, 0, 0, 0 }; 

uint64 pow[] = { 1/ 0 }; 

int na = 0; 

int i; 

cpy (a, _a, n) ; 
if (!a[l3 ScSc !a[2] && ia[3] ) 
{ a [03 /= SCALE; 

cpy (_a, a, n) ; 
return ; } 
if (is_negative (a, n) ) 

{ /* negative*/ 

negate (a, n) ; 
na = 1; } 
/* shift 'b r until it is greater than 'a'*/ 
for (i = 0; i < 128 && leq (b, a, n) ; i += 2) 
{ shiftup (b, n) ; 

shif tup (pow, 2) ; 
shiftup (b, n) ; 
shiftup (pow, 2) ; } 
if (is_zero (pow, 2) ) 

{ runtime$error (MSG_FL0AT0VF) ; /* overflow*/ } 

while (leq (scale, a, n) ) 



Appendix, Page 19 of 43 



/ * 

QIPLG Attorney Docket No:500.0028-20us 



{ /* while still a divisor*/ 

while (Iss (a, b, n) ) 
{ shiftdown (b, n) ; 

shiftdown (pow, 2) ; } 
if Ua[l] && la[2] ia[3] ) 
{ /* down to a 64-bit integer?*/ 

a [0] /= SCALE; 
add2 (r, a) ; 
break ; } 
sub (a, b ; n); /* subtract largest power*/ 
add2 (r, pow) ; /* increment result*/ } 

if <na == 1) 

{ negate (r, 2) ; } 

cpy, (^a, r, 2U} f f / ( f i / / / / f f f f f / / ; / / f / / f t ,,,,,,,,,,,,,,,,,, / 



////////// 

u REALIP GET INTEGRAL P ORTION OF A REAL ,,,,,,, 

, t , f , i , / / / / / / / / / / / / / 7 / / / / / / / 7 ' t / / / ? 7 i / / / / / / 7 / ( t > ///////////'///'/'' 

// Brief description: get the integral portion of a real removing the fractional part. 
// Expected: 

// in = address of input real 

// out = the address of a 64-bit real (void pointer) 

// Result: void : out = integral portion of in >..,,,,,,, , ,,,,,,, , 

t / / / 7 / / ///////////////////////'/'/''''' t f / / / / / / . ■ / / / / / / / / / * f f f t ' ' > f ' 
void routine realip (real * in, real * out) 
{ out->fp =0; /* no fractional part*/ 

, , ? u , t ; > , i ? :Ft?\v C,^^^ ,,,,,,,,,,;,;/,;//;/////// 

WHOLE INT GET INTEGRAL PORTION OF A REAL 



1 //////""/////////////////////// n i i i i i i i i i ' i < ' i t ' ■' ///////////'/'''''' 

h& If Brief description: get the integral portion of a real removing the fractional part. 
f"| // Expected: 

5f // in = address of input real 

|3 // out = the address of a 64 -bit real (void pointer) 

// Result: void : out = integral portion of in 

t i i iTiit / t tf > / / / / / / / / / / / t ft i > ' ' ' ' ' > ' ' ' 1 1 1 1 ' 1 ' 1 1 1 ' 1 1 ' 1 1 1 ' 1 1 ' ' 1 ' 1 { 1 1 1 1 
^Jj void routine whole_int (real *in, real *out) 

i J , f e ;»H P / (i / n '/ /° w ; / } / / / / / / / / / L < < / " <-L > > < < I* < < " > > < < < * > < 1 1 > 1 ><<><<<<< 1 



CVTLR CONVERT LONG TO A REAL 



*J! 1 7 7 , , i , , { , , , i i t i t i f i i i i i i i i t t f i i * * ' i i ' i ' * i ii 1 ii i ' i t t i / i i i i ' f > > * * 1 i * * 

p3 // Brief description: convert an integer to a 64 -bit real 

% ff Expected: 

Jl r . // a = input integer 

W // b = the address of a 64 -bit real 

W , / R ? s / U } t // } f ///////// 7/ a / //////////////////// ////// ////// / / / / / / / / / " ' 

J?* int routine cvtlr (int _a, real * b) 
{ real a = _a; 
*b = a; 

ff; , , f^ygSV] i !><>>>>>>> < > < > > > > < < < > < * > < > > > > 1 f > 1 f >' 1 1 1 1 ' 1 '* ' 

lM 1 a C ONVERT REAL TO A LONG ——- r — rr - rr - i 

i I i i i i i i i t i i i i t t i { i f / i i i i 7 / / 7 i t i i t i i f i i i i i f f t i f 



I f Brief description: convert a real to a long with rounding 
// Expected: 

// a = address of input 64 -bit real 

// b = address of output integer 

{', , rr?Y; ,},,,,,),,*,! rffftTr}?, ,,,,,, >> >> >>>> > - ,,<,>><><■>< >> 

int routine cvtrl (real *r, int *i) 

{ *i = (int) r->ip; // get integer part 

if (r->fp 1= 0 // Non-ZERO fraction? 

&& (r->fp >= SCALE 12 II and we have to round it 
| | r->fp <= - (SCALE/2) ) ) 

{ if (*i < 0) // round down if negative 

*i = *i - It 
else 

*i = *i + 1; // round UP if zero or positive 

} 

/ / f^/V/l /////////////////// / /// / /////////////// ^ ^ ^ ^ ^ / / ; 7 7 7 ; 7 7 7 1 

n GETINT GET AN INTE GER FROM A REAL ,,,,,,/,,//////; ' 

7 / ///////////////////////////////////// ii ii i ii 1 ' * < i * > * i * > * < ' < * > 1 1 ' 1 
II Brief description: Get an integer from a real 
// Expected: 

// n reference to a real 

// Result: integer value 
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\ 1 ' ' 1 1 * ' ' ' ' ' 1 ' i f * ' f * t * * i i / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 
mt routine getint (real & n) 
{ int i; 

cvtrl (&n, &i) ; 

return i ; } 

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
// I P 

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
// Brief description: Truncates a real number at the decimal point and 
// returns the integer portion. 
// Expected: in - passed real number 

// Result: out - the integer portion of a real number 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
void routine ip (int *out, real * in) 

{ *out = (int) in->ip; /* trunctate a int64 down to an int*/} // ip 

' / y f 1 ' f f f ' ' f f ' t f t t i f t f i i t / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / f / / i i i t i 



CMP64 COMPARE TWO 64-BIT INTEGERS 

/ / / / 1 1 f 1 1 ' f 1 * t f i t t t i t / / / / / 7 7 / / / / / / / 7 / / / i i i / / / / / / / / / / / / / / / / / / / / / / / / / 
// Brief description: compare two 64-bit integers 
// Expected: 

// a address of 64 -bit integer 

// b = address of 64 -bit integer 

// Result: 0 = equal 

// -1 = less than 

// 1 = equal 

l 1 ' ' ' ' ' ' > I t ' l I l / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ; / / / / / / / / / / 

int routine cmp64 (void *a, void *b) 
{ if (Mint64 *) a *(int64 *) b) 
return 0; 
if (*(int64 *) a < *(int64 *) b) 
return -1; 
fasj. return 1 ; } 

// f f ^ / f f f / f ///////// t //////////////////////,///,/,/////,/ f f f , j f 



DHEX64 DUMP A 64-BIT INTEGER IN HEX 

5Aa? /////////////////////////////////////// / 7 / // // / .'/////////////////// 

_ n II Brief description: Dump a 64 -bit integer in hex. This is a debug tool 

// f or 54 -bit arithmetic. 
'ill II Expected: 

III // buf = dump hex characters to this buffer 

]Z II _s = start address of 64 -bit quadword 

W // Result: void 

* f i f t f t t f f i f i i f t t t i / / / / / / / / / / / / / / / / / / / / / / / / / / / / , / / / / / / / / / / / / / / / / / / / 
Ki static void routine dhex64 (char *buf, void *_s) 
Q { unsigned int *s = (unsigned int *) _s; 
iTi sprintf (buf, »%08X%08X» / s [1] , s[0]);} 

W ' ft* ' ' ft * > f I > t ' t I ' t t i t I f f t I I I I t I I I t I I I I I I I f ! I i l > l ! t t { ! ! I I l ! I i f I f ( { / 

U 1 // DHEX96 DUMP A 96-BIT INTEGER IN HEX 

L; ''''//'/////////////////////////////////////// / t ////////////////// / 

|y // Brief description: Dump a 96-bit integer in hex. This is a debug tool 
Q. // for 96-bit arithmetic. 
J"! // Expected: 

*™ II buf = dump hex characters to this buffer 

// _s = start address of 96 -bit quadword 

// Result : void 

'/'//'///////////////////////////////////////////////////////////// 
void routine dhex96 (char *buf, void *_s) 
{ unsigned int *s = (unsigned int *) _s; 
dhex64 (buf, s + 1) ; 

sprintf (buf + 16, "%08X", s [0] ) ; } 

t t i i f i i / / ii f / f / / / / // / / // // / / // // // // // // // // // // // // // // // // // / / / / / 

h DHEX128 DUMP A 128-BIT INTEGER IN HEX 

i st / / / / / / / / / //////////////// //// ///////////////////,//,/,,,,,,/,, , 7 
// Brief description: Dump a 128-bit integer in hex. This is a debug tool 
// for 128-bit arithmetic. 
// Expected: 

// buf = dump hex characters to this buffer 

// _s = start address of 12 8 -bit quadword 

// Result: void 

////'////////////////////////////////////////////////////////////// 
void routine dhexl28 (char *buf, void *_s) 
{ unsigned int *s = (unsigned int *) _s; 

dhex64 (buf, s + 2) ; 

dhex64 (buf +16, s) ; } 

////////// ///////////////////////////////////////////////////////// 
u RTOD CONVERT REAL TO DOUBLE 

~ it / / / / / t /////////////////////////////////// / / 7~ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 



// Brief description: Convert a real to a double. 
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// Expected: 

// a address of real 

// b address of double 

// Result: void 

' ' '* ' '* ' ' ' ' J *' *' ' ' ' ' ' 1 ' ' ' * t ' i ' * * ' * t t t t i t i t i t t i t t i t t i i i t i t i t i i t t t t t t 
void routine rtod (real * a, double *b) 

{ double d = a->fp; 

d /= SCALE; 

d += a->ip; 
*b = d;} 

/ 1 1 / ' 1 ' 1 f 7 / ' / 7 ' ' f ' 1 L i i / / / / / / L_L i f t i i i i i ( i i t t f i i i i i i t t t i t t i i ; t ( t i t t f t 



„ GETINT GET A DOUBLE FROM A REAL 

1 ' ' ' 1 f ' ' ' ' ' 1 ' ' ' 1 < > ft f ft ////////// ////////////////////////////// 7-7-7 7 



// Brief description: Get a double from a real 
// Expected: 

// n address of a real 

// Result: double value 

''l'' 7 ''. ''''''''''''''''''Z////''////////////////////////////////// 

double routine getdouble (real * n) 
{ double d; 

rtod (n, &d) ; 

return d; } 

1 t ( l 1 1 l l 1 1 t 1 i 1 i 1 1 1 1 1 1 t j } t { ) / / / / / // // // // // // // // // // // // / / ///////// 
// D1V128 DIVIDE TWO 128-BIT INTEGERS — 

's ', 1 ' 1 ' 1 1 1 ' ' f 1 ' ' ( ' ' ' 1 ' ' 1 1 1 1 1 f 1 1 ' 1 7 1 1 1 ( 1 1 1 1 1 } 1 11 i> 1 11 ~ ~< 7 '< * f > t / 7 1 / 

// Brief description: Divide two 128-bit integers. This is used for intermediate 
// arithmetic. 
// Expected: 

// _dvend = address of real dividend 

// _dvsor = address of real divisor 

// quote = address of real quotent (division result) 

// rem = address of real remainder 

// Result: dividend/divisor = quotent, remainder 

' ' ' ' ' y 7 ' ' ' ' / ' 7 ' ' ' ' ' ' '''''''/>/'////////////////'//////////,,//,/// / 
void divl28 (void *_dvend, void *_dvsor, void *_quote, void *rem) 
{ uint64 dvend[] = { 0, 0 }; 

uint64 sdvsor[] = { 0, 0 }; 

uint64 dvsor[] = { 0, 0 }; 

uint64 quote [] = { 0, 0 }; 

uint64 pow[] = { 1, 0 }; 

int i; 

cpy (dvend, _dvend, 2) ; /* get local copy of dvend*/ 
cpy (dvsor, _^dvsor, 2); /* get local copy of dvsor*/ 
cpy (sdvsor, _dvsor x 2) ; /* get local copy of dvsor*/ 
if (is_zero (dvsor, 2)) 

{ /* no divisor ?*/ 

runt ime$ error (MSG_DIVBY0) ; /* divide by zero*/ } 
/* * shift divisor until it is >= to dividend 

* that is, whilst the divisor is less than the dividend.*/ 
for (i = 0; i < 128 && lss (dvsor, dvend, 2); i += 2) 

{ shi f tup (dvsor, 2) ; 

shiftup (pow, 2) ; 

shiftup (dvsor, 2); /* shiftup twice to speed up*/ 
shiftup (pow, 2) ; } 
if (is_zero (pow, 2)) 

{ runt ime$ error (MSG_FLOATOVF) ; /* overflow*/ } 

/* * Here we subtract the largest possible divisor value 

* from the divident each time adding power to the quotent. 

* We are trying to reduce the divident as quickly as we can.*/ 
while (leq (sdvsor, dvend, 2)) 

{ /* while still a divisor*/ 

while (lss (dvend, dvsor, 2)) 
{ /* follow the divident down in value*/ 

shiftdown (dvsor, 2) ; 
shiftdown (pow, 2); } 
if (! sdvsor [1] && ! dvend [1] ) 
{ /* down to just 64 -bit integers ?*/ 

uint64 n[] = { 0, 0 }; 
n[0] = dvend [0] / sdvsor [0] 7 
dvend [0] -= n[0] * sdvsor [0] ; /* remainder*/ 
add2 (quote, n) ; 
break; } 
sub2 (dvend, dvsor) ; /* subtract largest dvsor*/ 
add2 (quote, pow) 7 /* increment quotent*/ } 
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if (rem) 

cpy (rem, dvend, 2) ; /* remainder*/ 
cpy (_quote, quote, 2); /* return result*/} 

t 1 11 ' 1 ' *t f t t f i t f f / // / / / / / // // // / / / // / / / // // // // // // // / / // / / / / / / / / / / 



MUL128 MULTIPLY TWO 128-BIT INTEGERS 

///////////////////////////////////////////////////////>////>//>/// 
// Brief description: Multiply two 128-bit integers. This is used for intermediate 
// arithmetic. 
// Expected: 

// _a = address of 64 -bit integer array (_a * _b) 

// _b = address of 64 -bit integer array (_a * _b) 

// r address of 64 -bit integer array (result) 

// Result: void 

//////////////////////////////////////////////////////;//////////// 
void mull28 (void *_a, void *_b, void *_c) 
{ uint64 a[] = { 0, 0, 0, 0 }; 

uint64 b[] = { 0, 0, 0, 0 }; 

uint64 c[] = { 0, 0, 0, 0 }; 

int i; 

cpy (a, _a, 2) ; 
cpy (b, _b, 2) ; 

/* * make 'b' the smaller of the two*/ 
if (lss (a, b, 4)} 

{ int 64 temp [4] ; 

cpy (temp, a, 4) ; 

cpy (a, b, 4) ; 

cpy (b, temp, 4) ; } 
for (i = 0; i < 128 && ! is_zero (b, 4); i++) 

{ 

I* if (b[0] & 1) 

IS shiftdown (b, 4) ; /* shift down*/ 

^ shiftup (a, 4); /* <-- shift up*/ } 

cpy (_c, c, 2);} 

; : ~ ///////////////////////////////////////////////>/////////////////// 
^ n DIVR DIVIDE TWO 128-BIT REALS 

:.f| ///////////////////////////////////////////////'/////////////////// 

W // Brief description: Divide two 128-bit integers. This is used for intermediate 

!y // arithmetic. 

II Expected: 

// _dvend = address of real dividend 

* // _dvsor = address of real divisor 

13 // quote = address of real guotent (division result) 

= s j // rem = address of real remainder 

II Result: dividend/ divisor = quotent, remainder 
Ids /////////////////////////////////////////////////////////////////// 
L s void divr (real * dvend, real * dvsor, real * quote, real * rem) 
iy. { real q; 

Q r96di vide (dvend, dvsor, &q) ; 

if (rem) *rem = 0; 
1 * quote = q; } 

///////////////////////////////////////////////,/////////////////// 

n REAL::OPERATOR-() NEGATE A REAL 

/ / / / // // / / / / ////// t / // // // // // // // // t t / / / // // // ,/ // // / / // / / / // // // / 
// Brief description: negate the current real as in -real 
// Expected: The current real class 
// Result : real 

///////////////////////////////////////////////,/////////////////// 
real real : : operator - ( ) 
{ real r = *this; 
negate_real (&r) ; 
return r; } 

///////////////////////////////////////////////,/////////////////// 

a REAL: :OPERATOR+=() ADD A REAL 

/////////////////////////////////////////////////////////////////// 
II Brief description: add a real to ourselves as in real += real 
// Expected: The current real class 
// Result: real reference 

///////////////////////////////////////////////>/////////////////// 
real &real :: operator += (real r) 
{ /* real += real*/ 

addr (this, &r) ; /* add real*/ 

return *this;} 

/////////////////////////////////////////////// j ////////////////// / 

u REAL::OPERATOR++Q INCREMENT A REAL 

///////////////////////////////////////////////>/////////////////// 
// Brief description: increment a real by one as in real++ 
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// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////' iii i i / i / / / / // // // // // 
real &real :: operator ++ (int) 
{ /* real++*/ 

ip++ ; 

return *this ; } 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i i / / / / L I i i i i i i L * i L i f 1 / 

// REAL::OPERATOR-=() SUBTRACT A REAL 

////////////////////////////// //////////////// / //////////////////// 
// Brief description: subtract a real from ourselves as in real -= real 
// Expected: The current real class 
// Result: real reference 

) t , i i f f i I I i I I t I t t t I I I I l I / l I I I I I I / t l l 1 f I l f 1 I 1 l ! t I I t ' ( f i f * t i ' i ' i > 1 1 * > 

real &real :: operator -= (real r) 

{ /* real -= real*/ 

subr (this, &r) ; /* subtract two reals*/ 

return *this;} 

/ / / / / / / / / / i / i i i / / / i / / i i / / i / / / / / / / / / / / / / / / / / / / / / / / / / i i L * t * / ' f ' 1 1 1 1 1 1 

REAL::OPERATOR~0 DECREME NT A REAL ________ — 

//////////////////////// /////////////////////// / // // // // // // // // // / 
// Brief description: decrement a real by one as in real-- 
// Expected: The current real class 
// Result: real reference 

iiitfiiiittfJtittitttttittitttttitiftittittitfiittttittiiiitififiti 

real fcreal :: operator-- (int) 

{ /* real--*/ 

ip "" ; 

/ / f 6 /^ 1 / 1 / }________* i i i i i i f i i i i i i i i i i i i i i i * ' * * i i ' ' 1 1 1 f f f * 1 1 ! 1 ' ' 1 t 1 ' f 1 { 1 ' 1 1 ' 



n REAL::REAL() CONSTRUCT A REAL 

///////////////////////////////// ////////'//'//' /////////////////// 
// Brief description: construct a real and initialize it to zero 
// Expected: The current real class 

// Result: real reference ,,,,,,,,,,,,,,, 
/////////////////////////////////////////////////'//'''''' ' 11 1 f 1 1 ' ' 

real:: real 0 

{ ip = 0; 

/ / f P / / ___\ t i i i i i t i t i f i i i t i i i t // t / t ////// i t ///// t // t t i i i ' //ft/ LJ—L ii'i t 
f REAL::REAL() CONSTRUCT A RE AL FROM AN INTEGER 

f ! I f ( t t , I ! i I t I I / 1 i / l i I I I f I I I I ( I t 1 f I I I / I t I / I l f I I //// t // f / f f /////// I / 

II Brief description: construct a real from an integer 
// Expected: The current real class 
/ / Result * real reference 

////////////////////////////////////////////////////////////'/////' 
real : :real (int i) 

{ /* construct a real from an integer*/ 

ip = i; 

, i ~t ? '(^ t t i / i t / / i / t t / / / / / t i / / / i / i i t / i t i f j i i ' i * 1 1 1 1 1 i f 1 1 1 1 1 ' ' ' 1 ' 1 ' ! ' 1 
n REAL::REAL() CONSTRUCT A REAL FROM A 64-BIT INTEGER 



/,/ /////////////////J/////////////////' ' ' '''''' ' t I i l I I l l I I I I I I I I 1 1 1 

II Brief description: construct a real from a 64 -bit integer 
// Expected: The current real class 

// Result: real reference ,,,,,,,,,,,, 
ifffifiiiftiiitiiii///////////////////////////t////i/////i////////t 

real:: real (int 64 i) 

{ /* construct a real from a 64 -bit integer*/ 

ip = i; 

?p = 9 J ,,,,,,,,,, , / / / / i f / // // // // // // // // // // // 



///////////////////////////// i / i L ii if i i i ' i t / 
u REAL::REAL() CONSTRUCT A REAL FROM A DO UBLE 

////// / / / / //////////////// i / 7 i i f i i i / i i i i i i i i * ' < (////// i //////// i t f i 
II Brief description: construct a real from a double 
// Expected: The current real class 
// Result: real reference 

i i i f i i i i i i i i i i i i / / t i i i / i / i / / i i / / i i / / / i / / i i i i i i i i J ? 1 f J ^ f 
real:: real (double d) 
{ dtor (&d, this) ;} ^ 

f*,*?^? 11 } ,//,//////////////////,////;<<<"<<<<<<"<<<<<<<<> 
real:: real (int64 whole, int fract) 
{ ip = whole; 



/ / // d / z^ 8 / ? / } / ////////////////////////// / ////// / / ' ' / ' ' / ' ' ' ' ' ' ' ' ' / ' ' 



fp = fract; } 



////////////////// i /// / //////// / / / i 



//////////////////// i i i i i i i ' i 1 / 



„ REAL:: OPERATOR^) SUBTRACT A REAL FRO M ANOTHER REAL „ 

f / / j / / / / / / / / / / / / / / / / / / / / / / i / t i i f / / 7 / / i i / i i i / / / / / t / / i / / f / i / / t t / / i / / i 



II Brief description: subtract a real from another real 
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// Expected: The current real class 

// Result: real reference ,,,,,, 
i f t , t t , i / / / t t / i t / / / / / / / / / / / / / / / / / / / / / / / / / / / / * f t t / t t i / / / / / / / / / / / / / / / 

real real :: operator - (real i) 

{ /* real - real*/ 

real r = *this; 

/ / f j^SH / X t / I I t f i t t t I I f i i t I I I ! I I I i f i i t t i t t t ! t ! I I l I i t i t t i i i i t i i i i t i i 

a REAL::OPERATOR+() ADD A REAL TO ANO THER REAL 

//////////////////////////// /////////////////// / // // // // // // // // // / 
// Brief description: add a real to another real 
// Expected: The current real class 
// Result: real reference 

t i t i i i t i i i t i i t i i t t t i t t t i t t i t i t i i t t i i i t t i t i i i t t t t t t t t t t f t t t t i t t t i i i t 
real real :: operator + (real i) 
{ /* real + real*/ 

real r - *this; 

return r += i ; } 

/ / / // // // / / / // // // // // / / / / / / / 



f ! j i i i i t t I i t t l l f i i t i i t I i i i i t i t t t i i i t i i 



n REAL::OPERATOR<() COMPARE TWO RE ALS FOR LESS THAN 

1 ////////////////// 7 //////////////////////////// / ////////>///////''' 
// Brief description: compare two reals for less than 
// Expected: The current real class 

// Result: real reference ,,.,,,,,,,,>,,, 

int real :: operator < (real i) 

{ /* real < real*/ 

/ / f 6 *?™/ 1 / 8 / / (t ; h ; 1S / / ^1 ') i i t i i i i i i i t i t t f f i t t ' < > > > > > ' * > / / / / / / L ' ' ft ' ' i i ' ' ' ' ' I 
u REAL::OPERATOR<=() COMPARE TWO REALS FOR LESS THAN OR EQUAL 

1 ///////////////// 1 i i i i t i i t i i i i i i i i t f f t < t i i t ' ft ft /////////////////// 
fsk j! Brief description: compare two reals for less than or equal 
Q // Expected: The current real class 
22 // Result: real reference 

tJ y, ,//////////////////// i t t t i t t i i t t i i i t t t i t i i f i * * t * ////////// lit / / / / 

int real :: operator <= (real i) 
W { /* real <= real*/ 

f , , ffygVfg '*) 1 1 1 1 1 1 1 1 / L 1 1 1 1 1 1 1 1 1 tit 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

^3 1 „ REAL: :OPERATOR<-() COMPARE TWO REALS FOR GREATER THAN ___ 

flf 1 / / / , //////////////////////////////////// ///// it i / / > / ////////////// / 
15 // Brief description: compare two reals for greater than 
W. II Expected: The current real class 

* . i j i i s / u / V / / / r f *\ r / e f Tff/ //////////////////////////////////////"///"" 

f3 int real :: operator > (real i) 

jVl { /* real > real*/ 

r return Iss this) ; /* if 'i' is less then we are greater*/} 

p» z / / / / / / / / / / / / / / / / / / / / / / / / / / / i i i i i i i i i i i i i i i i i i ' f 1 I i i / L i / I i / L i ' L L i i 1 

pi 1 n REAL::OPERATOR>=() COMPARE TWO REALS FOR GREATER THAN OR EQUAL 

m /////////// / / / i i f i t / i i i i i f i t i t i t i i i i i i i j i i i i i t i i t i t t / i t i i i i i i i i i i i i 

Q // Brief description: compare two reals for greater than or equal 
PI // Expected: The current real class 

m , fr^i 1 1 ^s^rm /,,//,////////////////////////////////////"/" 

int real :: operator >= (real i) 

{ /* real >= real*/ 

if (lss (this, &i)) /* we are less*/ 

return 0; 

return 1; /* we must be >=*/} 

/ / 7 / I !)' I I t I I I I I I I I I I I I t I I I II fill I I I t I I t I I I I I I I I I I I I I I I I I I I i t i t i i f t 



REAL::OPERATOR%() THE MODULUS REMAINDER OF A REAL AND AN INTEGER 



//////////////////////////// / 7 //I II 11 I I f i I I f I I J ///////// I i / / I / / / / / f 

II Brief description: the modulus remainder of a real and an integer 
// Expected: The current real class 

{', , , R r, u }v , / , r ?^M e f e / r ?T? ,,///////////////////////////>///'/'/'"/""" 

real real :: operator % (int i) 
{ real a = (*this / (real) i) ; 

uint64 scale [] = { SCALE, 0 }; 

divl28 (this, (real *> scale, this, 0) ; 

real r = *this - (real) (a * (int64) i) ; 

/ / f e / t y r / n / r / / //////////////////////////////// ft fill ii ii if i i i i i / / / / / I t / / 



u REAL: :OPERATOR*() MUL TIPLY TWO REALS ,,,,,, 

t , , i i , , i , ! i i i t I 1 I ! , I i I i i i i i i i i i i i i t f i > ' f ' ' i ii ii i i t i i f i t I I i i i i i ii i i i 



II Brief description: multiply two reals 
// Expected: The current real class 

// f Result : f f real^ reference f f g f f t , , t f f t t t , , t , , , , f , , , , , t , t , , , , , f , , , t , , t , t i t 
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real real :: operator * (real m) 
{ real r; 

mulr (this, &m, &r) ,- 
return r; } 



///////// 



i i i i i f i i t i i i t i i t t i i i t i t t t i i i t i i t t i t t t i < t i t i i i i i i i i i i i i i i i i 



. REAL::OPERATOR&() A BIT AND OF A REAL AND AN INTEGER r ^_______ _ 

////////////////////////////// ///////////////// ' /////// f f f / t f / iff/ f 
II Brief description: a bit and of a real and an integer 
// Expected: The current real class 
// Result: real reference 

II I / III I I I / I I I / I I / I f f / I I f * * f * ' t t I I I I I I I I i I I I I I I ft ft / / / / / / / / / / ' / ' ' ' ' 

int real : : operator & ( int m) 

f return (int) (ip & m) ; /* bitwise AND on integer portion only*/} ,,,,,,,,, 
} / ill// / i ///////// / / t /// t i ////// i //// / 11 < i-L t-L > > LJ. L > ' L t ' 1 > ' ! 1 > < ' ' ' f 1 

„ REAL::OPERATO R/() DIVIDE TWO REALS _____ 

/////////// I t I I I I I I I I l I l I t t I I t I i l I I t f I I i I I I ! I I I I I / / / / / / / / / / / / / ' ' ' ' ' 

// Brief description: divide two reals 
// Expected: The current real class 

, / / , r f^/ r / e f e / r ?^ ,,/////////////////////////'/"/"/"""""" 

real real :: operator / (real m) 
{ real r; 

divr (this, &m, &r) ; 



return r; } 

/////////// 



//////////////////////////////////// // i ' i 11 1 / 1 1 ' 1 1 ' 1 1 ' 1 1 



III 



. R EAL::OPERATOR*=() MULTIPLY TWO REALS AND ASSIGN ,,,,,,,,,,, 

///////////////////////////////////// ///// ^ / / ' ' 7 t t /////// i / * 7-7 / / / / / 



/ l t t / l i / l l / l t i i I I I I I I I / t t ' t 1 



II Brief description: multiply two reals and assign 
// Expected: The current real class 
// Result: real reference 

( / I i i i i i i i i i i i i i t i i i i i t i i i i ' i i t f 1 i 1 * 1 1 1 
real real :: operator *= (real r) 

, { , f*/*?" 7^V/ ???? /V/* ^ / / / / / / / / LJJL1 > < < > * > LA > > ' " " < < IU < < /// ////// / / L < 

» REAL::OPERAT OR/=() DIVIDE TWO REALS AND ASSIGN ,,,,,,,,//, 

///////////// —i //////////////////////// /// ' ' ' ' 7 ' i / i /'/ f ' 
II Brief description: divide two reals and assign 
// Expected: The current real class 

, / / /^/^"Tf? ///////////////////////////'////'/"""""'" 

real real :: operator /= (real r) 

i , v^ff^fr, ???? /,') ! X , //////////////////////"/ " 

// CMPR^COMPARE TWO REALS FOR EQUALITY ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 
// Expected: The address of two reals 
// Result: -1 if a < b 
// 0 if a == b 

y / ////// /\ /f A ^ / ///////////////////////// ' ' ////// ' ' ' ' ' ' ' ' ' ' ' ' ' / ' 7 

int cmpr (real * a, real * b) 
{ if (*a == *b) 

return 0; 
if (*a < *b) 

return -1; 

return 1; J ^ ^ ; , / , ,,,,,,,,,,,,,,,,,,,//;///;////// / 



REAL::OPERATOR==() COMPARE TWO REALS FOR EQUALITY 



. , f ///////// /"/ / / / / / / / / 7 / / t / / 7 / / / / / 7 / / 7 / / 
// Expected: The current real class 
// Result: real reference 

t i t , i t i i I t I ! I I I I I I I i I l 1 I I i I i i i t > f i 1 1 * * 1 > 

int real: : operator -- (real r) { return ip == r.ip &_ fp == r 

/////// I I I I I I I I f I t I I I I I I I I i I I ! I I i ' I * 1 f 1 1 



f //////// t 



/ I I I / I I I I I 

fp;} 

I I I I I I I I I I 



t I / 

/ t I 
t I f 



//////// 



//////// 
I I I I I I I I 



u REAL::OPERATOR!=() COMPAR E TWO REALS FOR INEQUALITY 

i/iii/ii/ ii i / / /////// // 7 / / / / / /////// / 7 7 / 



// Expected: The current real class 
// Result: real reference 

/%///////// i i i i i i i i i / / i i i i i / / i i t i i i / i i i / 
int real: : operator != (real r) 
{ return I operator == (r);} 

//////////////////////// t /// /////// i i f 
REAL::OPERATOR!() LOGICAL NOT OF A REAL (IS ZERO) 



////////// 



////////// 



//////// 



//////// 



, t l I I I / I I I 1 I / i I l ! I I I t I / I I / I I I I l l t I * I * ( ' l 

II Expected: The current real class 
// Result: real reference 

/////// i // i i i i / i i i i i i i i i / i i i t / t i / / i i / i i i 
int real : : operator ! ( ) 

/ , f??? 1 / 1 / 1 ? " /' f / P / } / ///////// y /// y ///////// / / 



I I / I I I // I I 



I I I I I / I I / I 



//////// 



//////// 
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u REAL::OPERATOR»=() RIGH T BIT SHIFT OF A REAL AND ASSIGN ,,,,,,,, 

//////////;///////////////////////'/'' '/'''''''' ' ' ' ' ' ' ' ' ' ' ' 1 ' ' ' ' 1 ' ' 
II This function has not been implemented {it is not used) 
// Expected: The current real class 

, fTf?;, , fff^rm ,,,,,,,,,,,,,,,,,,,,,,,•..<:.<>>•">">>>">>>> 

real real :: operator (int) 

{ f return * this ;} / , f , t t f t/ /t f t /t t t f , t ,,,,,,, t , t f f f , it :f f/ f/ t ,,,,,,,,,,,, i 
u REAL::OPERATOR»=() LEFT BIT SHIF T OF A REAL AND ASSIGN ,,,,,, 

i / / 7 /////////"//////////////////////// > 1 1 > > * * > < ■ 1 ' J 1 ' 1 1 ' 1 1 1 1 1 ' 1 1 ' ' ' 

II This function has not been implemented (it is not used) 
// Expected: The current real class 

, frfft, , , r ??, r , e * e ?r?i ,,,,,,,,,,,,,,,,,,,,,,,,,,,•,.'.<<>>>>>><>>">< 

real real: : operator <<= (int) 

, { / X ^ r ?^f')f > U-UU-L t * * * * * > < > > ' * > > * > > /////// UL i ' ' UL UL > * ' < L ' > < 1 < 1 1 ' 1 1 



« REAL ::OPERATOR»=:() LEFT BIT SHIFT OF A REAL __ „ — — , 

////////////////////////////// //////// ' ' ' ' ' ' ' f ' ( ' 1 ' ' ' 1 ' 1 1 ' ' 1 1 1 1 1 



II This function has not been implemented (it is not used) 
// Expected: The current real class 

, , , , r f9, r , e ^J n ^ > . > >> >> ><>> > " " " ,>>>>>>><>>>><> ' 

real real :: operator « (int i) 
{ real r = *this « i; 

/ / ^TLull i , t i 1 1 1 i i * i i i * i >* i >** * ' ''*'*''< '''''''''''' ''''''' ''''''' * 



a REAL::OPERATOR»=() RIGHT BIT S HIFT OF A REAL 

/ / 77 / / / / / / / / / 7 //////////////// / / / / / * ' ' * > > < < ' ' ' * ' TTT-rn-, < i / / / / / / / / / / 
// This function has not been implemented (it is not used) 
a Expected: The current real class 

? , rvf?;, , r ?**?r?i ,,,,,,,,,,,,,,,<,,•,,,,,.<,<■<">"><<>>>""" 

15 real real :: operator » (int i) 
yg { real r = *this >> i; 
"Z return r; } 

// for rtoa() ascii digit table 
M #define TABLE_SCALE 10000 

#define TABLE_DIGITS 4 /* number of '0's*/ 

Btatic # char ntab [TABLES C ALE] [TABLE^DIGITS] ^ ,,,,,,,,,,,,,/,/////// 



> ' " i ^ A 7""-/ N / IT / n i /,//////////////// n T-r-r-rjy , > i n n > > n n n i i i , > ? / / / / / / / > 
f«s // Brief description: Initialize numeric table with ascii values 
H // The table is initialized with the least significant bytes 
W // first. For example, entry ntab[l] would be 

iyi // ntaMl] [0] = 'l'j // least significant digit first 

Q II ntab[l] El] = '0'; 

*w // ntabtl] [2] = '0'; 

C;| // ntab[l] E33 = '0'; 

m II and entry number ntab [1678] would be 

a ntab [1678] [0] = '1'; // least significant digit first 

// ntab [1678] [1] = '6'; 

// ntab [1678] [2] = »7' ; 

// ntab [1678] [3] = '8'; 

/"/ / / R ? s / u /V / rffi / / / /////////////////////'/"""' " " " ' " ' ' 1 " * * * * 1 1 * 

void routine math_init () 
{ char ascii_buf [323 ; 

char f ormat_st ring [32] ; 

int i, n; . 
sprintf (format_string, »%%0%dd», TABLE_DIGITS) ; /* build format string*/ 
for (i = 0; i < TABLE_SCALE ; i++) 

{ sprintf (ascii_buf, format_string # i) ; /* get ascii digits*/ 

for (n = 0; n < TABLE__DIGITS ; n++) /* move those digits into table*/ 

ntabEi] [n] = ascii^buf [TABLE_DIGITS - (n + 1)3; /* least significant first*/ 

, , , )] / / ,////////////////// ''/// / UL > > '-> L-L LJ. tU it * L > L I > t < 1 < ' t 1 ' < > 1 f 1 1 



GET ASCII F RACTIONAL_DIGITS ,,,,,,,,/,/////;//////> 

////"// 7 ////////////////// t it > ' f > > 1 ' 1 1 1 / 1 / ' f 1 ' ' ' 1 ' 7 ' ' ' 7 ' 1 1 1 



II Brief description: Get ascii fractional digits. 
// Expected: 

// a = address a character array for ascii digits 

// i = index into above character array 

// n = binary integer to convert 

// p = scale fractional digits 



Appendix, Page 27 of 43 



QIPLG Attorney Docket No:500.0028-20us 



// digits = number of precision digits 

// Result: i new index into a ,,,,,,,,, ,,,,,, , 

, /f/ 7 / / / / //////////////////////////////////'//'•' //////// / / / / / / / / / / / 
inline int routine get_ascii_f ractional_digits (char *a, int i, int n, 

int fract_digits, int digits) 

{ int table_digit, index; 

if (n < 0) /* calculate -1*/ 

n = -n; /* as 1*/ 

while (fract_digits > 0) 

{ if (n < TABLE_SCALE) 

{ /* already less no need for a divide*/ 

index = (int) n; /* use the number as the index*/ 
n - 0; /* no more number*/ } 

else 

{ index = n % TABLE__SCALE ; /* remainder*/ 

n /= TABLE_SCALE; /* new n*/ } 
for (table_digit = 0; table_digit < TABLE_DIGITS ; 
table_digit++, --fract_digits) 
if (fract_digits > 0 && fract_digits <= digits) 
a[i++] - ntab [index] [table_digit] ; } 

i , , , h , as " c " t srf±i ffigz >>< 

1 j ,^^™TWrf? > 

// Brief description: Get ascii whole integer digits. 
// Expected: 

// a = address a character array for ascii digits 

// i = index into above character array 

1^ // n binary 64-bit integer to convert 

□ I', , ^V; t??™, ^ f?? f, ,,,, ma / / / / / / ,///////////////////////// " " 

J- inline int routine get_ascii_whole_digits (char *a, int i, real * _n) 

W { int index, table_digit, skip_leading; 

ijj int64 n = _n->ip; /* integer portion*/ 

,f* do 

is f { if {n < TABLEJSCALE) 

{ /* already less no need for a divide*/ 

p| index = (int) n; /* use the number as the index*/ 

n = 0; /* no more number*/ } 

else 

si { index = (int) (n % TABLE_S CALE ) ; /* remainder*/ 

n /= TABLETS CALE ; /* divide by TABLE_S CALE * / } 
a[i++] « ntab [index] [0]; /* always at least 1 character*/ 
if (n) 

{ /* more to come*/ 

for (table_digit • 1; table_digit < TABLE_DIGITS ; 
table_digit++) 

p a[i++] = ntab [index] [table_digit] ; /* convert all digits*/ } 

?li else 

{ /* skip leading '0's on last divide*/ 

for (skip_leading = TABLE_DIGITS - 1; skip_leading ; 
--skip_leading) 
if (ntab [index] [skip_leading] != '0') 

{ for (table_digit = 1; table_digit <= skip__leadmg; 

table_digit++) 
a[i++] = ntab [index] [table_digit] ; 
break; 

} 

} } 

while (n) ; 

return i;) /* get ascii whole digits*/ 

return^ j f f f f / / / f f / / f y i f * f { f f f / / { f f / f f ,,,,, ,,,,,,,, f ,,,, i ,,,,,, t , / 



w 

5 s 



2 ~f 



tl RTOA REAL T O ASCII _ _ . . . . . ■ - ■ - . , . . . , , , , 

//////////////////////////////////// / ' ' ' ' 7 7 7 7 1 ' < ' i ' < i ' > < > < 1 > 1 1 1 f 1 f 1 

If Brief description: Convert a real for ASCII output. 

// This routine is almost identical to rtoa. . .except that it fills 

// in realinfo with information about the converted number. Informaticn 

// includes stuff like: whole_digits , f ract_digits , neg or not, ... 

// Expected: 

// in address of 64-bit real 

// _desc = address of output descriptor 

// f digits = number of fractional digits to display 

// flags = (2 = implied decimal point) ??fix this?? 

// Result: real_info is filled 
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, t i , , I } l , I f I i t I t f t I f ! ! I i l I I l l I I l I I I t I I I I I I I ' I ' 1 > 1 1 1 * { 1 ' ( 1 1 1 1 1 1 * 1 1 1 1 

void routine rtoa (real * in, void *_desc, int f digits, realinfo * rinfo, 

int flags} 
{ DESC_S *desc « (DESCJS *) _desc; 

char *r = desc->dsc$a_pointer; 

int len = desc->dsc$w_length; 

int skip_trailing, i = 0/ 

char buf [100] ; 

char *a = buf; 

real n; 

int f; 

int s = 0; 

int p; 

int is_implied_dec = (flags & 2) ; 

fill (sizeof (realinfo), rinfo, 0); /* zero all of the realinfo structure*/ 
n = *in; 

if (i senega tive_real (&n) ) 
{ negate_real (&n) ; 

s = 1; } 
rinfo ->neg = s; 

f = n.fp; /* get fractional part*/ 

/* * we build the string from left to right 

* in reverse order fraction. whole digits*/ 
if (fdigits | | f 1= 0) 

{ /* fraction first right to left*/ 

for (p = SCALEJDIGITS ; p < f digits; p++) 

a[i++] = '0'; /* fill with «0' digits first*/ 

i = get_ascii_fractional_digits (a, i, f, SCALE_DIGITS / fdigits); 
H~ for (skip_trailing =0; 

skip trailing < i && a [skip_t railing] == '0'; skip_t railing*--) ; /* skip »0's*/ 
g rinfo- >fractdigits = i - skip_t railing; /* number of fractional digits ignoring trailing '0's*/ 

O i f ( ! i s_imp 1 ied_de c ) 

yg a[i++] = '-'; /* add a decimal point*/ } 

;f| if {n.ip != 0) 

w f { /* get the whole ascii digits (all of them) */ 

tfl int new_i = get_ascii_whole_digits (a, i, &n) ; 

jfjj rinfo- >wholedigits = new_i - i; /* number of whole digits*/ 

i = new_i; } 
if (s) 
r« . if (i > 0) 

-f*\ a[i++] = '-'; 

rinfo->rlen = i; /* length of number (relative to 1)*/ 

yj a [i] = 1 \0 • ; // for easy debug 

Li* for (p = 0; p < len - i; p++) // blank fill result 

Hi r [p] = » 

for (; p < len; p++) 
fj { // fill result with answer 

si* r[p] = a[--i]; // from reverse order 

W }} // rtoa 

char* routine rtoa (real *n, int d) 
{ static char buf [28] ; 
DESC_S desc[l] ; 
realinfo info; 
desc->dsc$a_pointer = buf; 
desc->dsc$w_length = sizeof (buf) ; 
rtoa (n, desc, d, &info, 0) ; /* real to ascii*/ 

, , L L I > > > < < < > < > > < < > > < < ' > > >-L , i i i , t i * * i i * i * ' i * ' ' ' ' < ' ' < ' ' 11 11 11 

\u NAT CNVOUT CONVERT A REAL FOR ASCII OUTPUT 



m 



f ; rvrTT,:) Trrrrrrr?, r; < , , < > , * > m n - n m , , - - m t > - 

II Brief description: Convert a real for ASCII output - 
// Expected: 

// " val = address of 64 -bit real 

// result = address of ASCII buffer 

/ / digits = number of digits to display 

, ^VM^/m^TT)'/ ,,,,/,,//////////////////////////////'/"' 

int routine math_cnvout (real * val, void *result, int digits) 
{ DESC_S *desc ■ (DESC_S *) result; 

char *dest = desc->dsc$a_jpointer; 

int rlen = desc->dsc$w_length; 

char tbuf [48]; //?? why 48 

int len = 38; //?? why 38?? 
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€1 



y* 



int offset; 

DESC_S trap (tbuf , len) ; 
real info rinfo; 

fill (rlen, dest, 1 '); // blank-fill result 
rtoa (val, &tmp, digits, &rinfo, 0); 

offset = len - rlen; I 111 what if RESULT buffer is bigger than ours? 
if (offset < 0) 
offset = 0; 
copy (rlen, &tbuf [of f set] , dest); 
return TRUE;} 



///////////// 



i i i i i i i t i t ; i i t i i t i i / i i i i t i ' i 1 1 
/ / / / / / / / / / / / / / / / / / / / / / i / / / i i ' * 



iiit t / / / / iii i i i i f t i i i i i i 



« NEW_CNVOUT CONVERT A REAL FOR ASCII OUTPUT 

i i t i t i i i i i i t t f i t i i i t i i i i f i i i ' * i 1 1 1 ' ' 1 
II Brief description: Convert a real for ASCII output. 
// Expected: 

// val = address of 64-bit real 

// result = address of ASCII buffer 

// digits = number of digits to display 

// flags - (2 = implied decimal point) ??fix this to be a symbol?? 

// Result: realinfo structure is filled and returned 

y t f f , i i i i i t t I I I I I ! t I t l I I I t t I I t * I t t * 1 

static int routine new_cnvout (real * val, void *result, int digits, 

realinfo * rinfo, int flags, int max_len) 
{ DESCJS *desc = (DESC_S *) result; 
char *dest - desc->dsc$a_pointer; 
int rlen = desc->dsc$w_length; 
char tbuf [48]; //?? why 48 

int len = 38; //?? why 38?? 

int offset; 

DESC_S trap (tbuf, len); 

fill (rlen, dest, » '); // blank-fill result 
rtoa (val, fctmp, digits, rinfo, flags); 

offset = len - rlen; //?? what if RESULT buffer is bigger than ours? 
if (offset < 0) 
yj offset = 0; 

Jfi if (maxj. en < rinfo- >rlen) { fill(rlen, dest, ■*•); 

£1 return FALSE; 

ry } 

ffi copy (rlen, &tbuf [of f set] , dest) ; 

a * return TRUE;} // new_cnvout 

^a* const int cvtnumlen =32; 

H static int ten_power[] = { 1, /* 0*/ 

W 10, , /* w 

100, /* 2*/ 

1000, /* 3*/ 
10O00, /* 4*/ 

100000, /* 5*/ 

1000000, /* 6*/ 
10000000, /* 7*/ 

100000000, /* 8*/ 

loooooQOOo , , , , , , , , ft,,,,,,,,,,,,,,,,, ,,,, ,,*""""""""" 



n SCALE$REAL SCALE A REAL TO A NUMBER OF FRACTIONAL DIGITS 

i i , i , i i i , i i i i t j i t i i t i i t t * i tt 1 1 ( * < 1 ' ' i 1 ' 



i 



/ / / / / t / / / / / / / / / 7 i / / / / / / / 



/ l l l / l l l l l l I I I / I III I iiiiiiii i i < ' ' ' ' ' ; * 

// Brief description: Scale a real to a number of fractional digits 
// Expected: 

// out = address of 64 -bit real 

// in = address of 64 -bit real 

// f r ac = number of fractional digits (scale to) 

// Result: void 

// Example: 

// scale 87.36 by 1 = 87.4 
// scale 87.36 by 0 = 87 
// scale 87.36 by -1 = 90 

z 8 / 0 ? 1 /*/ 8 / 7 ) 3 / 6 /^/"? r, 1 , 00 , ,//////////,///////////""/ " " " 1 < * ' 11 1 f 1 ' ' ' ( 1 ' 

void old$scale$real (real * out, real * in, int frac) 
{ int64 scale [] = { SCALE, 0 }; 

int64 one[] = { 1, 0 }; 

int neg = 0; 

if (is_negative_real (in)) 
{ negate_real (in) ; 

neg =1; } 
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if (frac < -SCALE_DIGITS) 

frac = -SCALE_DIGITS; 
if (frac > SCALE_DIGITS) // our maximum precision 

frac = SCALE_DIGITS; 
/* * here we calculate the scale*/ 
if (frac < 0) 

{ /* -i == o.l*/ 

int64 pow[] = { ten_power [-f rac] , 0 }; 
mull28 (scale, pow, scale) ; } 

else 

{ int64 pow[] = { tenjpower [frac] , 0 }; 

divl28 (scale, pow, scale, 0); } 
if { I is_negative {scale, 2)) 

{ ~ " /* not negative*/ 

int64 n[] = { 0, 0 }; /* whole 'n'*/ 

int64 r[] = { 0, 0 }; /* remainder 'r'*/ 

int64 a[] = { 0, 0 }; 

int64 b[] = { 0, 0 }; 

a[0] = in->ip; 

scaleup (a) ; 

b[0] = in->fp; 

add2 (a, b) ; /* join ip & fraction*/ 

divl28 (a, scale, n, r) ; 

shiftdown (scale, 2) ; 

if (leq (scale, r, 2)) 
add2 (n, one) ; 

shiftup (scale, 2) ; 
ij* mull28 (n, scale, n) ; 

" n scale [0] = SCALE; 

W scale [1] = 0; 

P divl28 (n, scale, a, b) ; /* split the 128-bit integer*/ 

k f§ out->ip = a[0]; /* to integral part*/ 

*J out->fp = (int) b[0]; /* fractional part*/ } 

Msi if (neg) 

« W/v^ 1 , ,///,/// — ////// ,,,,,,, 



| 1 \ RpUNDSR^ROU^AREALTO A NUMBER OF FRACT;ONAL D.G.TS 

%a if Br ief description: Round a real to a number of fractional digits 

Si // Expected: 

// out = address of 64 -bit real 

H // in address of 64-bit real 

lij // frac = number of fractional digits (scale to) 

yL // Result: void 

hi // Example: 

W If Round 176.357143 to 2 fractional digits = 176.36 

Q // round 176.357143 to -1 fractional digits = 180 

m // scale 87.36 by 1 = 87.4 

// scale 87.36 by 0 = 87 

// scale 87.36 by -1 - 90 

// scale 87.36 by -2 = 100 

, , , ,////////////////////////'///''/'' / / / / ' " ' / " ' ' ' 1 1 1 ' 1 1 1 ' 1 1 1 1 ' 7 
void routine round$real (real *out f real *_in, int frac) 
{ real rin = *_in; 

real *in = &rin; 

int scale = SCALE; 

int rem; /* remainder*/ 

int neg - 0; 

if (is_negative_real (in)) 
{ negate_real (in) ; 

neg =1; } 
if (frac < -SCALE_DIGITS) 

frac = -SCALE_DIGITS; 
if (frac > SCALE_DIGITS) /* our maximum precision*/ 

frac = SCALE^DIGITS ; 
/* * positive r fraC then only scale the »fp' part*/ 
if (frac == 0) 

{ out->ip = in->ip; /* copy over ip part*/ 

if ( (SCALE » 1) <- in->fp) /* 'fp' part greater than a half*/ 

out->ip += 1; /* carry one*/ 

out->fp =0; /* no decimal digits*/ } 
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else if (frac > 0) 

{ out->ip = in->ip; /* copy over ip part*/ 

scale /= ten_j>ower [frac] ; 
rem = in->fp % scale; /* get remainder*/ 
out->fp - in->fp / scale; /* scale down*/ 
if {(scale >> 1) <= rem) /* is there a carry?*/ 

out->fp +=1; /* carry one*/ 

out->fp *= scale; /* and up again (to get rid of lower digits)*/ 

if (out->fp >= SCALE) 
{ /* overflow ?*/ 

OUt->fp -= SCALE; 
out->ip++; /* carry one*/ 
} } 

else 

{ /* scale the 'ip' part*/ 

int64 ip - in->ip; 

if ( (SCALE >> 1) <- in->fp) /* 'fp' part greater than a half*/ 

ip += 1; /* carry one*/ 

out->fp =0; /* zero result fp part*/ 

scale = tenjpower [-frac] ; /* a ten_power*/ 
rem = (int) (ip % scale); /* get remainder*/ 
out->ip = ip / scale; /* scale down*/ 
if {(scale » 1) <= rem) /* is there a carry?*/ 

out->ip += 1; /* carry one*/ 

out->ip *= scale; /* and up again (to get rid of lower digits)*/ } 

if (neg) 

/ / / / ^f/VfM W/l / / / / /////// / / / / / UL < < >,< i <> i <<>>>> <> > t-L < < " < * > * LI 
n SCALE$REAL SC ALE A REAL TO A NUMBER OF FRACTIONAL DIGITS ——^ 

If Brief description: Scale a real to a number of fractional digits 
// Expected: 

// out = address of 64 -bit real 

// in = address of 64 -bit real 

// frac = number of fractional digits (scale to) 

// / Result: / void^ / / / / / t f f / / f 9 9 f / / § i f f f f f f t 9 f f f f f t , ( f , , , , , , , , , , , , , , , , , , , 

void scale$real (real * out, real * in, int frac) 
{ round$real ((real *) out, (real *) in, -frac);} 
char *dotptr; 
char *cvtnumptr; 

r^v^ zin^ff^h l £±i¥h ? e > r <°< d > lSl &}f^h * ! > ,,,,,,,,,,iiiii<" 

n FL TO STRING FRAC r~n~i~T~r~n~r~7~i~i 



, 7 } l l , / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' i i } i i i < 1 f 1 1 //,///// 
II Brief description: Convert a float to a string 
// Expected: val - value to be converted to a string 
// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of fractional digits 

i f , / , / / / / z^ 1 * 9 ?/ 

int fl_string_output (int *len, char *outptr) ; 

int routine new_f l_string_output (int *len, char *outptr, int flags); 

int routine f l_to_string_f rac (int *len, char *outptr, real * val, int fdigits, 

int flags) 
{ char tmpbuf [cvtnumlen] ; 

DESC_S tmpdsc[l3; 

real rvalue ; 

int status; 

real info rinfo; 

int is_implied_dec ■ (flags & 2) ; 
int *s - (int*) &val->ip; 
*outptr = '\0' ; 
if (len[0] <= 1) 
return 0; 

set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 
pad = (flags & cvt_pad) != 0; 
rvalue = *val; 
s = (int*) &val->ip; 

new_cnvout (&rvalue, tmpdsc, fdigits, &rinfo, flags, *len) ; 
if (tmpbuf [0] '*') 
return 0; 

zeros = fdigits - rinfo. fractdigits; 
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dotptr = fctmpbuf [cvtnumlen - rinfo.rlen] ; // point us to the first data 
neg = rinfo.neg; 
if (neg) 

dotptr = &dotptr[l]; // skip sign (??dme?? fix this later) 
sci = FALSE; 

leader = rinf o.wholedigits; 
period = leadzero = 0; 

frac = f digits; // rinfo.f ractdigits; 

//?? fix & 2 to be " IMPLIED_DECIMAL " 
if (is_implied_dec) 

status - new_f l_string_output (len, outptr, flags) ; 
else 

status = f l_string_output (len, outptr) ; 
if (! status) 

fill (len[0], outptr, '*•); 

return 1;} // f l_to_string_f rac 

////////////////////////////////////////////////////////;'///////// 

// CVT OUT D F 

>>>//////////>/////////>>////////////////////// / /////////////////// 
II Brief description: Convert real to floating ASCII output. 
// Expected: 

// fval = address of 64 -bit real 

// dlen = fill length 

// desc = pointer to output descriptor 

// Result: status from fl to string_frac 

////////////////// ~i rt ////////////J////////////////'////////''''''' 
int routine cvt_out_d_f (real * fval, int dlen, void *desc) 
{ int status = SUCCESS; 

int len; 

char *text; 

char numbuf [100] ; 

split$desc (desc, &len, &text) ; 

round$real (fval, fval, 0) ; 

status = f l_to_string_f rac (&len, numbuf, fval, 0, 0) ; 
if (! status) 

return MSG_CVTERROR ; 
// len - = 1; 
fill (dlen, text, > '); 
copy (len, numbuf, text + dlen - len) ; 

return status;} // cvt out d f 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i t f f t t * ' * t t ' / / / 1 

y REAL CVT OUT D F 

///////////////////////// ////////////////////// / / // // // // // // // // // 
// Brief description: Convert real to floating ASCII output. 
// Expected: 

// fval = address of 64 -bit real 

// dlen = fill length 

// desc = pointer to output descriptor 

// Result: status from fl to string frac 

///////////////// / /~//~////~7//////////////////////////////'//''//'' / ' 
int routine real_cvt_out_d_f (real * fval, int dlen, void *desc) 
{ int status = SUCCESS; 

int len; 

char *text; 

char numbuf [100] ; 

char *p = numbuf; 

split$desc (desc, &len, &text) ; 

round$real (fval, fval, 0); 

status = f l_to_string_frac (&len, numbuf, fval, 0, 0) ; 
if ('status) 

return MSG_CVTERROR ; 
fill (dlen, text, 1 ' ); 
if (dlen < len) 

{ p += len - dlen; 

len = dlen; } 
copy (len, p, text + dlen - len) ; 
return status ; } 

/////////// ///////// j /////////// / / / // / ///// ^ f / t t ////// ////////// f / / 
NAT CVT STR FLT CONVERT ASCII FLOATING STRING TO REAL 



/ / / /~/ / 7 1 i 7 ////////////////// //////////////////////// / / ' ' t ' t ' 7 ' ' / f ' 1 
II Brief description: Convert ascii floating string to real 
// Expected: 

// result = address of result 64-bit real 

// dlen = length of ascii text 
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t f , , / ,///,////////////////////////' / ' ' " " " / / " " it < <" < i 

int routine math_cvt_str_f It (real * result, int len, char *text) 
{ real n = 0; 
int neg = 0; 
int e = 0; 

for \i = 0; i < len && textli] && text [i] <= ' »; i++) ; /* skip leading*/ 

if (text[i] == •-■) 

{ neg = 1; 

i ++ ; } / 

for (; i < len && text [i] text [i] <= ' »; ; /* skip leading*/ 

if (i >= len) 

{ *result = n; 

return 1; } 

if (! (tolower (text[i]) >= '0' && tolower (textli]) <= '9')) 
{ // not numeric 

switch (tolower (text [i] ) ) 
{ case 1 e ' : 

case 'x' : 
case ' + 1 : 
case 1 . 1 : 
breaks- 
default : 

return 0; 
} ) 

/* * get the integral part first*/ 

for (; i < len && text [i] >= '0' && text [i] <= '9'; i++) 
{ n.ip *= 10; 

n.ip += textti] & OxOf; } 
if (i < len && tolower (textti]) == 'e') 
{ e = 0/ 

for (i++; i < len && text [i] >= '0' && textti] <= '9'; 
{ e *= 10; 

\S e += text[i] & OxOf; } 

while (e--) 

n.ip *= 10; } 
if (i < len && text [i] == '.') 
8 { int f = 0; 

2*1 int p = ++i; 

if (i < len && (textli] < '0' || text [i 3 > '9')) 

return 0; 
for (; 

i < len && i < p + SCALE_DIGITS && text [i] >= '0' 
SW && text [i] <= »9' ; 

{ f *= 10; 

f += textti] & OxOf; } 
/* do not round ascii input*/ 

/* if (i == p+6 && textti] >= '5' ScSc textti] <= '9')*/ 
/* f += 1; */ 
for (; i < p + SCALE_DIGITS ; i++) 

f *= 10; 
n.fp = f; } 
if (i < len && tolower (text[i3) == 'e') 
{ e = 0; 

for (i++; i < len && textti] >= '0' && text[i] <= '9'; i++) 
{ e *= 10; 

e += textti] & OxOf; } 
while (e--) 
{ n.fp *= 10; 

n.ip *= 10; 
if (n.fp >= SCALE) 

{ n.fp -= SCALE; 

n.ip++? 

} 

} } 



s7s 



m 



if (is negative_real (&n) ) 

{ runtime$error (MSG_FLOATOVF) ; /* overflow*/ } 



if (neg) 

negate_real (&n) ; 
* result = n; 
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return 1 ; } 
int routine ator(real *a, char *b) 
{ / ptm^th^Btomja, ^trlentb^b) t\ 

„ NAT CVT STR_FLTJMPDEC CONVERT ASCII FLOATING STRING TO REAL 

////"//"// 7 ///////////////// / / t i / / / / / t < / ' ' ' > < ' ' 1 ' /////'/''' / / / / 
// Brief description: Convert ascii floating string to real 

// 

// Expected: 

// result = address of result 64-bit real 

// dlen « length of ascii text 

// text = ascii input text 

/ / ?W/ /////////////////////" " " " " " : ' ' ' ' ' ' " ' 1 1 1 1 1 1 ' 1 1 1 

int routine math_cvt_str_f lt_impdec (real * result, int len, char *text # int f digits) 
{ real n = 0; 
int neg = 0; 

int i ; a. / 

for (i = 0; i < len && text [i] && text [i] <= ' '; ; /* skip leading*/ 

if (text[i] =- ' - 1 ) 
{ neg = 1; 

i++ ; } , 

for (; i < len && text [i] && text [i] <= ' '; i++) ; /* skip leading*/ 

if (i >= len) 

{ *result - n; 

return 1; } 
/* * get the integral part first*/ 

for (; i < len-f digits && text [i 3 >= '0' && text [i] <= '9'; i++) 
{ n.ip *= 10; 

n.ip += text[i] & OxOf; } 
int f = 0; 



fl int p = i; 

^ if (i < len && (text[i] < '0' || text [i] > '9')) 
w return 0; 

•v0 for (; i < len && i < p + SCALE_DIGITS && text [i] >= '0' 

% && text [i] <= '9' ; i++) 

ff { f *- 10; 

f\J f +« text[i] & OxOf; } 

111 for (; i < p + SCALE _DIGITS ; i++) 

^ f *= 10; 

* n.fp = f; 

l'3 if (is_negative_real (&n) ) 

s*j { runtime$error { MSG_FLOATOVF ) ; /* overflow*/ J 

ff if (neg) 

negate__real (&n) ; 

HI * result = n; 



return 1; 



//////////// / ' ' / / ' ' ' ' / / ' / ' ' ' 7 ' 1 f 1 > 1 ' 1 1 1 1 1 1 1 



Bl 1 j DTOR CONVERT A DOUBLE TO REAL / 

l« 1 / 7 7 / / / / / / 1 / / / / / / / / 7 i / / / t i / t t / / / / / / 7 7 / / / / ' / / / ' 



/ / f/ffttlfitt titiiit''''' ; • 

// Brief description: Convert a double to a 128-bit real 
// Expected: 

// a * address of input double . 

// b address of output 128-bit real 

// f Result: / void / / f / / y y / y y / / ; / / f ; / / , , , , , , , , , , , , , , , , , , , , , ifi , , ,,,,,,/,/ / 
void 

routine dtor (double *a, real * b) 
{ double d = *a; 

double my_d; 

int 64 ip; 

int fp; 

if (d > LARGE S T__NUMBER | | d < -LARGEST_NUMBER) 

runtime$error ( MS G_FLOATOVF ) ; 
// handle rounding now by adding in 5*10^-8 
if (d < 0) 

my_d = d - 0.000000005; 
else 

my_d = d + 0.000 0000 05; 

ip = (int 64) my_d; 

my d -= ip; 

fp~= (int) (my_d * SCALE) ; // this is the same code I last sent you 

fp = (fp/10)*10; // truncates at SCALE-1 digits 
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if (ip ScSc fp < 0) /* if ip then*/ 

fp - -fp; /* fp always positive*/ 

b->ip = ip; 
b->fp = fp; 

/ / / } / //////////////////////// > ' ' ' L > * 11 1 * / / ' ' 1 11 > 1 ' ' ' ' f ' 1 ' ' ' 1 ' ' 1 ' - - ' 1 

u MULR MU LTIPLY TWO REALS ,,,,,,,//,;/,;///;//// ' 

/,/////////////////////////'/''' ~ ' ' ' ' ! 1 / 1 I l 1 1 J ! j. ' 
// Brief description: Multiply two 128 -bit integers. This is used for intermediate 

// arithmetic. 
// Expected: 

// _a = address of real {_a * _b) 

// _b = address of real (_a * _b) 

// r = address of real (result) 

// f Result: f void / /f /fti tf f ft f f ,,,,,,,,,,,, ,//,/,,//.,//////////////'// ' 
void tnulr (real * _a, real * _b, real * _c) 

\ , f^r^f^M' r?% T C ) ! X !_!_/_/_/ ,,,,,,,,,,,,,,,.■,,,>•>>><>"">"" 

I j CyT$QUAD$FLT CONVERT A 64-BIT QUADWORD TO REAL 

// Brief description: Convert a 64-bit quadword to a real. 
// Expected: 

/ / result = address of real 

// p_quad = address of 64 -bit quadword 

// / Result : / statiw f f ,,,,,, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,,/,,/,/// / / / / / / 
int routine cvt$quad$flt (real * result, int *p_quad) 
{ result->ip = (*(int64 *) p_quad) ; 

result->fp =0; // Clear fractional portion ++dme++ 

/ / fy??/ 1 /! 1 lu > < < > ' < < > * < < > < * * > < * * f 11 * ' 1 1 1 1 1 1 ' 1 1 1 1 1 1 1 1 1 — I 

FL TO STRING — J 



; TY,, t , , , , \ > / 7 i > i , s / / / / / / / * * > > / / / / " " rn > > n - 7-7 n m n n n < > 



i i i i i i i 



7 / I I iiiftflfiii*'* tt,( '' . 

II Brief description: Convert a float to a string 
// Expected: val - value to be converted to a string 
d // " outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of significant digits 

M / / / / / / /W / ///////////////////"""'"/' ' ' .V " ' ' ' ' ' ' ' ' ' 1 1 f 1 ' 1 1 1 ' 

& int routine fl_to_string (int *len, char *outptr, real * val, int digits, 
3 int flags) 

f 1 { char tmpbuf [cvtnumlen] ; 
v " DESC_S tmpdscEU; 

W int status; 

1^ real info rinfo; 

ks pad = (flags & cvt_pad) 1= 0; 

!.:f set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 

Q new_cnvout (val, tmpdsc, digits, &rinfo, 0, *len) ; 

m if (tmpbuf [0] '*■) 

return 0 ; 

zeros ■ digits - rinf o. f ractdigits; 

dotptr = Stmpbuf [cvtnumlen - rinfo. rlen] ; // point us to the first data 
neg = rinfo. neg; 

if < ne 9) , . -, _ x 

dotptr = &dotptrtl]; // skip sign (??dme?? fix this later) 

sci - FALSE; 

leader = rinf o. wholedigits ; 
period = leadzero = 0; 
frac = rinfo. f ractdigits; 
cvtnumptr = tmpbuf; 

status = f l_string_output (len, outptr) ; 

retur V tatus;} / ,_uu d i &Ell5±E< ,,,,,,,,,,,< <<><><> ">>">> " 



: / F , L / ?T?' < W ® /°"/ T / p /" T / ; 



// Brief description: Convert a float to a string 
// The result looks like: [- 3 nnnnnnn. f f f f f f f 
// Where [-3 is either a space or ' - 1 
// nnnn = whole number, f f f f =f raction 
// (notice how the decimal point is included) 
// Expected: val - value to be converted to a string 
// outptr - an address of conversion. 

// len - maximum length of conversion 
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// digits - the number of significant digits 

' t ' f / / / / / / / f W/ 

int routine f l_string_output (int *len f char *outptr) 
{#define $save_lit (ch) \ 
{ \ 

if ( (outptr - outbegin) > maxlen) \ 

return 0; \ 
outptr [0] = ch; \ 
outptr++ ; \ 

} 

#undef $save 

#define $save(tlen, txt) \ 
{ \ 

if ( (outptr - outbegin) > maxlen) \ 

return 0; \ 
copy(tlen f txt, outptr); \ 
outptr += tlen; \ 
} 

char * outbegin; 
int maxlen; 

memset (outptr, 0, *len) ; /* null fill result*/ 

maxlen = len[0] ; 

outbegin = outptr; 

if (neg) 

$save_lit ('-') 
else 

if (pad) 

$save_lit C ' ) ; 
O if (leadzero) 

ft* $save_lit ('0'); 

W $save (leader, Scdotptr [0] ) ; 

'JI if (frac == 0 && leader == 0 I leadzero) 

>M $save_lit CO'); 

% if (frac > 0) 

*M { $save_lit ('.'); 

HI $save (frac, &dotptr[l + leader]); } 

ffi if (pad) 

^ $save_lit ( ' * ) ; 

* *len = outptr - outbegin ; 

Q if (maxlen < *len) 

s' i return 0; 

U ^is^vj ' u ^i^m 



ify.v.vr?, ? 



/ / / 7/ / i i f j i f i t t i i t i t > ' * ' ' • • • ' • \ ■ 

// Brief description: Convert a float to a string 
T*i // The result looks like: [- 3 nnnnnnnf f f f f f f 
jflj // Where [- ] is either a space or ' - ■ 
// nnnn = whole number, f f f f =f raction 

// (notice how the decimal point is implied but not included) 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len - maximum length of conversion 

/ / digits - the number of significant digits 

,",/////, , £ W 8 / w 7 W^m/^ 1 /^/ /,////,//////////////////"'"""" 

int routine new_f l_string_output (int *len, char *outptr, int flags) 
{ char * outbegin; 
int maxlen; 

int is_implied__dec - FALSE; 
if (flags & 2) 

is_implied_dec = TRUE; 
maxlen = len[0] ; 
outbegin = outptr; 
if (neg) 

$save_lit ('-') 

else 
if (pad) 

$save_lit C '); 
if (leadzero) 

$save_lit CO*); 
$save (leader, fcdotptr [0] ) ; 
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if (frac == 0 && leader == 0 && Ileadzero) 
$save_lit CO'); 

if {frac > 0) nn 
{ if (is_implied_dec) //??dme?? as WHy no » ; " on line below?? 

$save (frac, fcdotptr [leader] ) // just the decimal digits... nc ».« 

$save S !frac + l, fcdotptr [leader] ) ; // fract+1 to include the } 

if (pad) 

$save_lit ( 1 ' ) ; 
*len = outptr - outbegin; 
if (maxlen < *len + is_implied_dec) 

return 0; 

return 1 • ) // new f 1 string output 



, F /ST? / / / A; / > i * * /////////////// n , ? n i i 7 i :— n , , > > 

II Brief description: Convert a float to a string 
// Expected: val - value to be converted to a string 
// " outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of significant digits 

/"/////// / f w g ? / /////////////////////""''"/'/' ' ' ' ' ' /////"""""' 
int routine real_f l_to_string (int *len, char *outptr, real * val, int digits, 

int flags) 
{ char tmpbuf [cvtnumlen] ; 
DESCJS tmpdsc [1] ; 
int d = 0; 
L.i int i; 

m* pad = (flags & cvtjpad) 1=0; 

set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 
p math_cnvout (val, tmpdsc, digits); 

-.fS if (pad) 

^ { *outptr++ = ' • ; 

m d - ii } 

for (i = 0; tmpbuf [i] ; d++) 
3«1 *outptr++ = tmpbuf [i] ; 

?W if (pad) 

09 { 

*outptr++ « 1 ' ; } 
*len « d; 
e = 0 ; 

lij frac = 0; /* leader - zeros;*/ 

s '*" sci = 0 ; 

P leader = 0; 

fy period - leadzero = 0; 

0 , , £S/Vj , ,,,,,,,,,,,,,,,,,,,,,,,,,,, 

W 1 ; AS T? z 8 /^/'"/ 0 / f ^7 vTT, 7 7 / 7 n , / ; / / / / / » 7 , , , 71 , , , Tl , , Tl Tl n n , , , , , , » I » 
// Brief description: Convert a float to a string 
// Expected: val - value to be converted to a string 
// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of fractional digits 

// flags: 

// if ( ! fl_to_jstring_frac (&numlen, buffer, (int*) &number, 
II form->form$f raclen, 0)) 

/ , / H*?™ / /////// / / / / / / ' ///////// " > << <<<<<<<<" >< " " << > " 11 1 

int routine real_f l_to_string_f rac (int *len, char *outptr, real * val, 

int digits, int flags) 

{ char tmpbuf [cvtnumlen] ; 
DESC_S tmpdsc [1] ; 
char *e = outptr + *len; 
char *p = outptr; 

char fill = '0'; /* default*/ 

int i, n; 

if (len[0] <= 1) 

return 0; 
pad = (flags & cvt_jpad) I- 0; 
set$desc (tmpdsc, sizeof (tmpbuf), tmpbuf ) ; 
math cnvout (val, tmpdsc, digits); 
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if (pad) 

*p++ = 1 ' ; 

for (i = 0; tmpbuf [i] && p < e && tmpbuf [i] 1» * .'; 

*p++ = tmpbuf [i] ; 
if (digits) 

{ if (tmpbuf [i] ! = '.') 

*p++ = ' . ' ; 
else 

*p++ - tmpbuf [i++3 ; 
for (n = 0; n < digits tmpbuf [i] && p < e; i++, n++) 

*p++ - tmpbuf [i] ; 
for (; n < digits && p < e; n++) 

*p++ ' 0 ■ ; 
fill = 1 •; } 

else 

{ if (tmpbuf [i] 1= 1 . M 

{ 

fill = ' } 
for (; tmpbuf [i] && p < e; i++) 
*p++ = tmpbuf [i] ; } 

if (pad) 

*p = ' 
if (p < e) 

{ /* move to end of string*/ 

while (p > outptr) 

*--e = *--p; 
while (e > outptr) 

U *" e - fil1 '* ) 

0 , . . fSVzj J± ^^"^ , , t u-u-l 

o \ 'fy*™**??7™? t >><" 

ip // Brief description: Convert a float to a string 

% II Expected: val - value to be converted to a string 

W // outptr - an address of conversion. 

y3 // len - maximum length of conversion 

ft\ // digits - the number of significant digits 

IS Z 7 , / / / / / / /"W 1 / 

int routine real_f l_string_output (int *len, char *outptr) 
{ char * out beg in; 
int maxlen; 
maxlen = len[0] ; 
outbegin = outptr; 
if (neg) 
fy $save_lit C-M 

y if (pad) 

fl| $save_lit ( ■ » ) ; 

if (leadzero) 

$save_lit COM; 
$save (leader, &dotptr[l]); 
if (frac == 0 && leader == 0 && \ leadzero) 
$save_lit CO 1 ) 
else 
if (frac > 0) 

{ $save_lit CM; 

period = 0; 
if Csci) 

while (e < 0) 
{ e++ ; 

$save_lit CO')} 
$save (frac, &dotptr[l + leader]); } 
else if (Isci) 

{ while (e - leader > 1) 

{ e--; 
$save_lit COM} } 
if (period) 

$save_lit CM; 
if (sci) 

$save (4, &cvtnumptr [cvtnumlen - 4] ) ; 
if (pad) 



3! 



3 I 
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$save_lit ( ' » ) ; 
*len = outptr - outbegin; 

return 1;} // f l_string_output 

#include "tmath.h" 

/////////////////////////////////////////////////////////////////// 

u NAT INT SIN 

///////////y//////y///////////y////y////y////y///////////////////// 
II Brief description: Returns the SINE of an angle specified in radians 
// Expected: in - passed real an angle 
// Result: out - sine of given angle 

////////////////////y///yy///////////////////////////////////////// 
double mth$dsin {double *) ; 

coroutine int math_sin (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *} in); 
b = mth$dsin (&a) ; 
*(real *) out = b; 

return 1;J // math_sin 

///////////////////////////////////////////////>/////////////////// 

n NAT INT COS 

/ / / // >/ // />/ / i / / // // / / / > t /////////////////// i //,////////////// i i t i i 
II Brief description; Returns a cosine of an angle. 
// Expected: in - passed real angle 
// Result: out - cosine ofd an angle. 

///////y///////////////////////////////////////. '/////////////////// 

double mth$dcos (double *) ; 

coroutine int math_cos (double *in y double *out) 
{ double a; 
real b; 

a - getdouble { (real *) in) / 
b = mth$dcos (&a) 
*(real *) out = b; 

return 1;} // math_cos 

/////////////// / / / / / i //////// / / / / // // / / / // / / / / / / / / / i / i / t / / / i / it i ( / t 

| u NAT INT TAN 

/////////////////////////////////////////////////////////////////// 
// Brief description: Returns a tanget of an angle. 
// Expected; in - passed real angle 
// Result: out - tanget of an angle 

/////////////////////////////////////////////////////////////////// 
double mth$dtan (double *) ; 

coroutine void math_tan (real * in, real * out) 
{ double a, b; 
real r; 

a - getdouble (in) ; 
b = mth$dtan (&a) ; 
dtor (&b, &r) ; 

*out = r;} // math tan 

I t l l I I I I l I t l I i t I t l I l l I i I t ! I I I I J I I I J i l I f I I f f ( f t i / / t t f / i t i t / i t i i i i i i t 

n NAT INT ASIN 

///////////////////////////////////////////'//// / //////////// / / / / t / i 

II Brief description: ASIN{x) returns the angle whose SIN is x. 
// Expected: in - passed real value 

// Result: out - angle 

/////////////////////////////////////////////////////////////////// 
double mth$dasin {double *) ; 

coroutine void matn_asin (double *in, double *out) 
{ double a; 
real b; 

a - getdouble {(real *) in); 
b = mth$dasin (&a) ; 

*(real *) out = b;} // math_asin 

//////////////////////y///y/y/y// / / // // // // / // / // // // / / / / / / // / / / / / / 

h NAT INT ACQS 

/iff/// / / // // // // /y//////y/////// / / / / // / / / / / / // /////////////////// / 
II Brief description: ACOS (x) returns the angle whose COS is x 
// Expected: in - passed real value 
// Result: out - angle . 

/////////////////////////////////////////////////////////////////// 
double Tnth$dacos (double *) ; 

coroutine void math_acos (double *in, double *out) 
{ double a; 
real b; 

a ■ getdouble ({real *) in) ; 
b = mth$dacos (&a) ; 

*(real *) out = b;} // math_acos 

i i i t t i i i i i i i i t i i i / i i / i i t ( i i i i i i i i t t f i t t i i i f i i t f t t / i i i i i i i i i t i f i i i ' ' 
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7 / N A T / 1 N / T / / A 7/ A / N / ////////////// > mrrrrnTTT-n n / / n n * * / / n n / / / / / / / / / 

// Brief description: ATAN(x) returns the angle whose TANGENT is x. 
// Expected: in - passed real value. 

, 777?! • , 7 7^}V, /,/////////////////'//'"//'/""""""""'" 

double mth$datan (double *) ; 

c_routine void math_atan (double *in, double *out) 
{ double a; 
real b; 

a = getdouble {{real *) in); 
b = mth$datan (&a) ? 

, , : ( ^»^ 7*>\ IL ,,,, t <">•> 



' W, ~ n , , n , 7 , i 7 / n i , i , , n n n n , m , , n m , , , m , m , n 

II Brief description: SINH(X) returns hyperbolic sine of X. 
// Expected: in - passed real value. 

, rr?Y;,,°^, -,, ^pj^h f^j^rrr^r, ,,,,,,,,,,,.,,,,,,,,,<■,,>■>.<■ 

double mth$dsinh (double *) ; 

c_routine void math_sinh (double *in, double *out) 
{ double a; 
real b; 

a - getdouble ((real *) in); 
b = mth$dsinh (&a) ; 

, , as r, b ;l , , , , , , m a t N s ^ h , , ,,,, ,,,,. ,,,,,,,,,,..,....>> 

' 7^7 W, f° S " 7~7~7~7~7~7~7~7~7~r7~7~7~7~r*7~7~7~7~7~7~7~7~7~7~7~7~i 

f I Brief dexcription: 

// COSH{X) returns hyperbolic cosine of X. 
// Expected: in - passed real value. 

, ,rjv- r, ^^,7/ a /? a /■r^™v , ^/ ,,,,,,,,,,,,,,,,,,,,,,,,,,,, 

double mth$dcosh (double *) ; 

c_routine void math_cosh (double *in, double *out) 
{ double a; 
real b; 

m a = getdouble ((real *) in); 

ffi b - mth$dcosh (&a) ; 

r ^tW,^ IL 



yj // Brief description: TANH (X) returns hyperbolic tangent of X. 
// Expected: in - passed real value. 

, / / Tti ' ?, ?*T, d /?* u f ,/,////,////////////"/"""' 

double mth$dtanh (double *) ; 
W coroutine void math_tanh (double *in, double *out) 
|U { double a; 

real b; 

a = getdouble {(real *) in); 
b = mth$dtanh (&a) ; 

. . «g r, b A I 1 , , m ath ^ n h L-Lj—i—L 

U^A 1 , ///////// r> r-m n n n r~, n > n , > > i , , , i < < / r , , > / / / / / / / ' ' > > > ' ' > 

II Bri4f description: // Expected: value - numeric expression seed = global seed value 

, 77??! Tf*r?, rfrfrrrrr, r?,*?, fr?*?, !',,,,,,, , v r7r?, rrrri 1 ?*,-, 

void realip (real *, real *) ; 
real mth$ random (int *) ; 
void math_rnd (real * value) 
{ real f = mth$ random (&seed) ; 

, , r?}™.-,v),, ,,,,,, Z£J 5SzgS - 



II Brief description: Returns the sign of a number 
// Expected: in - passed numeric expression 
// Result: out - sign of a a numeric expression 
// = +1 if expression > 0 

// = -1 if expression < 0 

ii,,,, A j £ , TfTf^, ,,,,,,,,,,,.,.,,.•.><"> " 

c_routine int mth$sgn (int *) ; 
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coroutine int realsgn (void *) ; 
void math_sgn (real * in, int *out) 
{ real n = 0; 

t out [01 # - $cmp$real On, /& n) ,-} // math.sgn ^ ,,,,,,,,,,,/,,/ 



* NAT INT ABS ^ — — ,,,/,////;/////;/////// 



jp:Ks 



Sis \ 



II t\ r\ I | |i | /A U w 

f { I I t ! , , I I I I i i I I f 1 1 I I I I I I ' I * I * * t I 1 1 * ' 1 ' ' 1 ' 1 ' 

ff Brief expression: 
// Returns the absolute value of a numeric expression. 
// Expected: in - passed real expression 

, f-r^; ,T, 7 ^ f/ 8 ? 1 "/ 6 / ^™ ° f / 

double mth$dabs (double *) ; 
void math_abs (real * in, real * out) 
{ double a; 
real b; 

a = getdouble ((real *) in); 
b = mth$dabs (&a) ; 

t , : w/;^ t II r^^i , , , , — 



7 r£ri ?f i i-mm rrrrm rri-rrrrri n , n n ,,/,,;, ,,,,,, 

// Brief description: Returnes the larger of two reals X and Y. 
// Expected: x - real value 
// y - real value 

, ^// ffr, , 1 7TF, < 7/ t r? f^' 7/ / ,/////////////-/////////////////" 

void math_real_max (real * x, real * y, real * out) 
{ if ($cmp$real (x, y) < 0) 
*out = *y; 
else 

,,,,:<??,-,?',> f ±^jS±Sj-LJu ' ' 



■; , N , A , T , ^^VA ,,,,,,,,,.,>,,,, >>>>"'4-"""' "" 



// Brief description: Returns the lesser of the two reals X and Y 

// Ejected x P - passed real value // y - passed real value // Result: out - lesser of two 



reals X and Y / / / / / / / / / / / / / / # g f f f f / , , f , , , , , , , , , , , , , , , , , , , , , , , , , , ,,,,,// / 
void math_realjnin (real * x, real * y, real * out) 
m { if ($cmp$real (x, y) < 0) 

*out = *x; 
* else 

,,,,:°ff,-,7',> '^^T^Vft 1 , , , , ,,, ,,,, 



" 1 j ^ f ?* rTTTTTTTTT, 7 , , , V 7"7~ rT^ , //^ ' " V ' ' 'J ' ' 'l number ? 



I- ',/ description^ ReturnVa ^'^ia' re^ value // Expected: in - passed real number 

ry // / fT^v^s ,,,,,,,,,,,**<,<>*!<><<<<<><<<<><<<<"<>""" 

G double mth$dsqrt (double *) ; 
f_ void math_sqr (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in) ; 
b = mth$dsqrt (&a) ; 

r° u * r, b A, , , , , , , , L / 1 ■ 11 / ,,, ,,,,,,,,//// ^ ^/ " 



///////////////// / ^ ' / / / 7 



// BAaf descAption; ^tuins't^a natural 'logarithm of a real — // Expected: in - passed real value 

, / R ? s / U /V / / ° u ^/ ? a / t " a ' }°?? x ^> /',,,,,,,,">>>><><>>>>>"">""""""" 

double mth$dlog (double *) ; 
void math_log (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in) ; 
b = mth$dlog (&a) ; 

, , r, b /), , "> t , ,,,,,,,,,,.,,,,,...<>>> — 

l : / N , A // 7???, ,,,,,,,,,,,, 777,, ,,,,,,,,,,,,,,,,,,,.>•>,•.•>>•>>> 

// Brief description: Returns the base 2 logarithm of given real value 

,", fWW . ft , - 7 p ^ rT} 7f , , 77??} , r?, ; 2 < ,,,,,,,, 

double mth$dlog2 (double *) ; 

void math_log2 (real * in, real * out) 
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{ double a; 
real b; 

a = getdouble (in) ; 
b ■ mth$dlog2 (&a) ; 

, . zM > >> t > > > ' ' " ' ' ' '' ?ftL& 

j / N / A / T ; m / L /°/ G / 1 /°/ n-r-m > , > , , , i , t > i > • i > > > ttttt-} n n , n ? > > i n / / / / 7 / / / / / / 
// Brief description: Returns the common logarithm of the real value 

// Ejected: t in / - passed / real val^ // / Result S/ out - common logarithm / ,,,,,,,,,,,,,,,,, 

double mth$dloglO (double *) ; 
void math_loglO (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in) ; 
b = mth$dloglO (&a) ; 

*out = b;} // math_loglO 
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